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ABSTRACT 

We present a numerical simulation of the dynamical collapse of a nonrotating, 
magnetic molecular cloud core and follow the core's evolution through the 
formation of a central point mass and its subsequent growth to a 1 M & protostar. 
The epoch of point-mass formation (PMF) is investigated by a self-consistent 
extension of previously presented models of core formation and contraction in 
axisymmetric, self-gravitating, isothermal, magnetically supported interstellar 
molecular clouds. Prior to PMF, the core is dynamically contracting and is not 
well approximated by a quasistatic equilibrium model. Ambipolar diffusion, 
which plays a key role in the early evolution of the core, is unimportant during 
the dynamical pre-PMF collapse phase. However, the appearance of a central 
mass, through its effect on the gravitational field in the inner core regions, leads 
to a "revitalization" of ambipolar diffusion in the weakly ionized gas surrounding 
the central protostar. This process is so efficient that it leads to a decoupling of 
the field from the matter and results in an outward-propagating hydromagnetic 
C-type shock. The existence of an ambipolar diffusion-mediated shock of this 
type was predicted by Li & McKee (1996), and we find that the basic shock 
structure given by their analytic model is well reproduced by our more accurate 
numerical results. Our calculation also demonstrates that ambipolar diffusion, 
rather than Ohmic diffusivity operating in the innermost core region, is the 
main field decoupling mechanism responsible for driving the shock after PMF. 

The passage of the shock leads to a substantial redistribution, by ambipolar 
diffusion but possibly also by magnetic interchange, of the mass contained 
within the magnetic flux tubes in the inner core. In particular, ambipolar 
diffusion reduces the flux initially threading a collapsing ~ 1 M core by a factor 
^ 10 3 by the time this mass accumulates within the inner radius (~ 7.3 AU) 
of our computational grid. This reduction, which occurs primarily during 
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the post-PMF phase of the collapse, represents a significant step towards the 
resolution of the protostellar magnetic flux problem. 

Our calculations indicate that a 1 M Q protostar forms in ~ 1.5 x 10 5 yr for 
typical cloud parameters. The mass accretion rate is time dependent, in part 
because of the C-shock that decelerates the infalling matter as it propagates 
outward: the accretion rate rises to ~ 9.4 M Myr" 1 early on and decreases to 
~ 5.6 M Myr -1 by the time a solar-mass protostar is formed. The infalling 
gas disk surrounding the protostar has a mass ~ 1CT 2 M at radii r ^ 500 AU. 
A distinguishing prediction of our model is that the rapid ambipolar diffusion 
after the formation of a protostar should give rise to large 1 km s _1 ), and 
potentially measurable, ion-neutral drift speeds on scales r <; 200 AU. 

The main features of our simulation, including the C-shock formation after 
PMF, are captured by a similarity solution that incorporates the effects of 
ambipolar diffusion (Contopoulos, Ciolek, & Konigl 1997). 

Subject headings: accretion, accretion disks — diffusion — ISM: 
clouds — ISM: magnetic fields — MHD — stars: formation — stars: 
pre-main-sequence 



1. Introduction 

It is generally accepted that most of the star-formation activity in our galaxy takes place 
through the gravitational collapse of molecular cloud cores (e.g., Mouschovias 1987; Shu, 
Adams, & Lizano 1987). It is, furthermore, believed that interstellar magnetic fields play a 
central role in this process in that their stresses are the dominant agent that acts against 
gravity to prevent, or delay, cloud collapse (e.g., Mouschovias 1978; McKee et al. 1993). 
This is embodied in the concept of a critical mass M cr i t , which in general takes account of 
both the magnetic and the thermal pressure contributions to the support of the cloud, but 
which, in the case that magnetic stresses dominate, reduces to M crit « O.lS^s/G 1 / 2 , where 
G is the gravitational constant and <3>£ is the magnetic flux that threads the cloud. Clouds 
whose mass M exceeds M crit are "supercritical" : they collapse on the free-fall timescale. In 
contrast, "subcritical" clouds (M < M crit ) can avoid collapse on the much longer ambipolar 
diffusion timescale. In the latter case, the neutrals gradually contract by diffusing inward 
through the ions and field, leaving behind a magnetically supported envelope and eventually 
forming a supercritical core that undergoes dynamical collapse (e.g., Mouschovias 1996). 

Because of the complexity of the problem — it involves solving the full dynamical 
equations of a magnetized, multicomponent (neutrals, ions, electrons, as well as charged and 
neutral grains) fluid that evolves over many decades in size and density in a nonspherically 
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symmetric manner (because of the presence of ordered magnetic fields and likely also 
rotation) — much of the progress in this area has been accomplished through the use of 
numerical simulations. One of the main efforts to simulate core formation and contraction 
due to ambipolar diffusion in magnetically supported molecular clouds has been carried 
out by Mouschovias and coworkers (e.g., Fiedler & Mouschovias 1992, 1993; Ciolek & 
Mouschovias 1993, 1994, 1995, hereafter CM93, CM94, CM95, respectively; Basu & 
Mouschovias 1994, 1995a, b). These studies followed the evolution of the core over six 
decades in density up to central densities ~ 3 x 10 9 cm -3 , where the assumption of 
isothermality starts to break down because of radiative trapping (e.g., Gaustad 1963; 
Hayashi 1966). This assumption had been adopted in the interest of simplicity: sophisticated 
and computationally intensive numerical techniques are generally needed to calculate the 
thermal structure of the gas during the opaque phase of protostellar evolution (e.g., Larson 
1969, 1972; Winkler & Newman 1980; Stahler, Shu, & Taam 1981; Boss 1984; Myhill 
& Kaula 1992; Myhill & Boss 1993). As a result of this restriction, the aforementioned 
calculations did not follow the collapse of the core to the time where a point mass - a 
protostar - is formed at the center, although they did obtain valuable information on the 
conditions leading to this critical event. In particular, by the time these simulations were 
terminated, the inner region of the core was collapsing dynamically and was characterized 
by neutral infall speeds ~ C (the isothermal speed of sound) and inward accelerations 
^ 0.3|g r | [where g r (r) is the local gravitational acceleration]. Furthermore, the thermal 
pressure, while remaining relatively unimportant in the envelope, came to exceed the 
magnetic pressure near the center. Basu (1997) derived a time-dependent, semianalytic 
solution that extended these ambipolar diffusion models up to the instant of point-mass 
formation (hereafter referred to as PMF []). He found that ambipolar diffusion continues 
to gradually erode the retarding magnetic forces in the inner core, making the collapse 
increasingly more dynamical (and the thermal-to-magnetic pressure ratio in the inner core 
progressively larger) as PMF is approached. 

The diminution of magnetic forces in the innermost regions of a collapsing core just 
prior to PMF suggests that one could gain some insight into the protostar formation process 
from previous studies of PMF in nonmagnetic, spherically symmetric, isothermal clouds. 
Analytic similarity solutions have uncovered two limiting behaviors: Penston (1969) and 
Larson (1969) found a solution in which, just before PMF, the infall speed approaches 
~ 3.3 C at all radii while the density scales with radius as r -2 , resulting in a spatially 



1 Creation of a central point mass was commonly referred to in previous papers as "core formation." In 
this paper we use the term "point-mass formation" so as not to confuse this process with the formation of 
an extended, magnetically supercritical, molecular cloud core. 
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uniform mass inflow rate ~ 29 C 3 /G (where G is the gravitational constant). Hunter (1977) 
extended this solution past PMF and showed that, immediately after the central mass is 
formed, the accretion rate onto the protostar increases to ~ 47 C 3 /G. In the other limit, 
Shu (1977) obtained a solution that is static prior to PMF (with the density distribution of 
a singular isothermal sphere, which also scales as r~ 2 ) and that takes on an expansion-wave 
character (with a constant mass accretion rate ~ 0.98 C 3 /G onto the protostar) following 
PMF. Hunter (1977) and Whitworth & Summers (1985) demonstrated that there are, in 
fact, infinitely many similarity solutions that span the range between the L arson- Penst on 
and Shu results, with the nature of any given solution being determined by the initial 
configuration of the cloud and the conditions at its boundary. Numerical simulations 
carried out by Hunter (1977) and by Foster & Chevalier (1993) confirmed the dependence 
on the initial and boundary conditions. In particular, it was found that the behavior of the 
central regions of clouds that are initially marginally stable to collapse approximates that of 
the Larson-Penston solution at the PMF epoch, although it was determined that the mass 
accretion rate onto the protostar declines at later times. It was, however, also found that 
the post-PMF evolution of clouds that initially have more extended envelopes approximates 
that of the Shu solution at late times. Since the initial and boundary conditions of real 
clouds are expected to depend on the detailed configuration and evolution of the embedded 
magnetic field, it is clear that one needs to incorporate magnetic field effects into the 
collapse calculations to adequately model the formation of protostars. 

There have been several recent attempts to calculate PMF following the collapse of 
magnetic interstellar clouds. Although they have all contributed to our understanding 
of the processes involved, their results were hampered by the adopted assumptions or 
approximations. For example, Tomisaka (1996) modeled clouds that had equal thermal 
and magnetic energy densities, so that they were not primarily supported by magnetic 
fields. This means that his model clouds were magnetically supercritical. This assumption 
is at variance with H I and OH Zeeman measurements of magnetic field strengths in 
molecular clouds, which are consistent with models of magnetically subcritical clouds 
(Crutcher et al. 1993, 1994, 1996). Li & Shu (1997) modeled PMF in self-similar, 
magnetic cores. They assumed that cores immediately before PMF can be represented 
by hydrostatic configurations of singular isothermal disks and that the magnetic flux is 
frozen into the neutrals. These assumptions are inconsistent with the above-cited results of 
numerical simulations and semianalytic solutions of the collapse of magnetically supported 
molecular clouds that undergo ambipolar diffusion (as well as with the simulations of 
thermally supported spherical clouds), which have found that the inner core regions collapse 
dynamically as PMF is approached (see also §3.2). As we show below, ambipolar diffusion, 
which plays a key role in bringing about the dynamical collapse, is generally important also 
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after PMF. Safier, McKee, & Stahler (1997) did include ambipolar diffusion in the modeling 
of magnetic cloud collapse and post-PMF evolution. However, their model is strictly 
spherical, and they neglected the effect of thermal pressure gradients in comparison with 
magnetic stresses. Furthermore, their formulation is quasistatic and does not incorporate 
the magnetic induction equation for the time evolution of the magnetic field. Li (1998) 
extended the Safier et al. model by including thermal-pressure and time- dependent terms 
and by adding the induction equation. This enabled him to follow the time evolution of his 
model cores (for r > 150 AU) even during the dynamical phases of the collapse. However, 
by retaining the spherical-symmetry assumption of Safier et al., his calculations were also 
unable to yield the geometry of the magnetic field lines. f\ 

The importance of ambipolar diffusion after PMF can be inferred by comparing the 
ambipolar diffusion timescale r AD = r/v D (where v D is the ion-neutral drift speed) and the 
gravitational contraction (~ free-fall) timescale r gr = (r/|g r |) 1//2 before and after PMF. In 
axisymmetric geometry, r AD /r gr ~ (r gr /r ni ) hb 2 in the inner flux tubes of a supercritical core 
(e.g., Mouschovias 1991), where //e(r) is the total mass-to-flux ratio at radius r (in units of 
the critical value for gravitational collapse) and 



the average collisional rate between ions of mass m; and neutral H 2 molecules of mass 
m H2 (~ 1.7 x 10~ 9 cm 3 s _1 for collisions between neutrals and Mg + or HCO + ions; 
McDaniel & Mason 1973); the factor 1.4 accounts for a 20% helium abundance by number. 
CM94 and CM95 found that, to first order for the late pre-PMF evolution of cores in 
disk-like clouds, the magnetic field B « 3 (r /r) mG (where r = 40 AU), rii « 0.1 cm -3 
(valid for neutral densities n Q ^ 10 7 cm -3 ; see Figs. 2b and 46 in CM94), the vertical 
column density cr n ~ 5(ro/r) g cm -2 , and the total mass M(r) k 6 x 10~ 3 (r/ro)M Q . 
For a disk-like cloud, the critical mass-to-flux ratio (M/$£)d,crit = (47r 2 G) -1 / 2 , where 
G is the gravitational constant (Nakano & Nakamura 1978). Therefore, one finds 
/ie ~ 2.7, \g r \ « GM(r)/r 2 k 2 x 10~ 6 (r /r) cm s~ 2 , and r ni « 7 x 10 9 s, which 
yields r gr « 2 x 10 10 (r/r ) s and r AD /r gI . « 20 (r/r ). It follows that r AD /r gr > 1 for 
r ^> 2 AU. Hence, as has already been demonstrated in the past, ambipolar diffusion 



2 The same is true for the spherically symmetric self-similar model of a collapsing magnetic cloud devised 
by Chiueh & Chou (1994), which, however, does not include ambipolar diffusion. It should be noted that all 
models that assume spherical symmetry do not satisfy the solenoidal condition V • B = on the magnetic 




(1) 



is the neutral-ion collision time. In equation ([!]), rt\ is the ion density and (crw) n i is 



field. 
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is ineffective as a dynamically collapsing core approaches PMF, and for r > 2 AU the 
magnetic flux can be considered frozen into the neutrals. (For reasons discussed in §2.2, 
we do not consider r <^ 5 AU.) Turning now to the post-PMF epoch, when the central 
point mass comes to dominate the gravitational field in the innermost flux tubes, we 
use the relation r AD /r gr w (1 — |a n |/|(7 r |) _1 (r gr /r ni ), where a n is the inward acceleration 



of the neutrals and g r is the total gravitational acceleration (see eqs. [fL4j] and |15| in 
§3.3). In this case, normalizing r as before, \g r \ ~ 4 x lO~ 4 (M cent /M Q )(r /r) 2 cm s~ 2 . 
Substituting again r ni rs 7 x 10 9 s (corresponding to n\ ~ 0.1 cm -3 ) , the above ratio 
becomes r AD /r gr « 0.2 (r/r o ) 3 / 2 (M /M cent ) 1 / 2 (1 - \a n \/\g r \y\ For \a n \/\g r \ in the range 0.2 
- 0.9, which corresponds to the period after PMF when the collapse becomes progressively 
more dynamical, one infers T AD /r gr ps (0.2 — 2)(r/r o ) 3 / 2 (M /M cent ) 1 < /2 . Hence r AD ^ r gr for 
-^cent ^ O.1M and r ^ 50 AU. [A similar result is obtained if one continues to use the 
pre-PMF relations and simply substitutes M cent for M(r).] This estimate indicates that 
decoupling of the gas and magnetic field by ambipolar diffusion should occur in the inner 
core regions after PMF. The physical reason for this is that the strength of the gravitational 
field in the weakly ionized gas near the origin is greatly enhanced by the appearance and 
growth of a central point mass, causing the neutrals to fall in more rapidly while the 
plasma and magnetic field are left behind. (The same basic reason — the appearance of a 
progressively growing free-fall zone around the origin following PMF — is also the cause 
of the increase in the mass inflow rate into the center at that epoch first discovered in the 
above- referenced nonmagnetic collapse calculations.) The foregoing conclusion is verified by 
a detailed calculation in § 3.3, where we show that ambipolar diffusion after PMF is, in fact, 
so efficient that it effectively decouples the neutrals and magnetic field in the innermost 
core region, with dramatic consequences for the subsequent dynamical evolution of the core. 
An alternative, yet equivalent, analysis of the effectiveness of ambipolar diffusion following 
PMF, based on the scaling of magnetic forces (particularly the magnetic tension force) after 
PMF, is given in Appendix C. 

In this paper we present a detailed numerical simulation of point-mass formation in 
a nonrotating, magnetic, dynamically collapsing protostellar core, properly accounting for 
the effect of ambipolar diffusion. Unlike earlier studies, we use an initial state that is 
consistent with the realistic models of core formation and collapse, as presented earlier 
by Mouschovias and coworkers. As we noted above, the simulations carried out by that 
group were terminated at densities where the isothermality assumption started to become 
invalid because of radiative trapping. Although a proper treatment of radiative trapping 
is indispensable for a complete treatment of the star-formation process, one can adopt 
a simpler approach that circumvents this difficulty by removing the region of radiative 
trapping from the active computation mesh. This is justified by the fact that the region 
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of radiative trapping (typically a few AU) is several orders of magnitude smaller than the 
characteristic core size (r corc « 0.1 pc). This region can therefore be considered to be 
effectively point-like, and one can proceed to calculate the formation of a central point mass 
within a gravitationally collapsing core and its effect on the subsequent evolution of the 
core by retaining the isothermality assumption. In adopting this approach, we note that the 
assumption of isothermality was also employed in previous studies of protostar formation in 
nonmagnetic cloud cores as well as in the more recent attempts to model PMF in magnetic 
clouds. 

The plan of the paper is as follows. In § 2 we review the main characteristics of the 
pre-PMF core-evolution calculations of CM93, CM94, and CM95 and outline the model 
modifications that we have implemented to extend the simulations beyond PMF. In § 3 we 
present the results of our calculations. We consider a typical model and follow the evolution 
of all physical quantities of interest through the PMF epoch. We also describe the formation 
(due to ambipolar diffusion in the innermost core after PMF) of a C-type shock and consider 
its propagation through the collapsing core. The appearance of such a shock as a result of 
field-matter decoupling was first pointed out by Li & McKee (1996), who proposed that 
the relevant field decoupling mechanism was Ohmic dissipation in the innermost regions 
of the core. Our study, however, reveals that ambipolar diffusion occurring outside the 
region of Ohmic dissipation is the main cause of magnetic flux decoupling after PMF. In 
§ 4 we discuss the structure of the ambipolar diffusion-mediated shock and show that our 
numerical calculations qualitatively reproduce the predictions of the simplified analytic 
model constructed by Li & McKee (1996). We present a quantitative comparison with their 
results and also address the issue of interchange instability in the post-shock region, first 
raised in their work, in light of our detailed computations. In that section we also discuss 
the observational implications of our simulation and briefly comment on the magnetic flux 
problem during star formation. Our results are summarized in § 5. 

2. Characteristics of Model Clouds 

2.1. Pre-PMF Phase 

For the pre-PMF phase of core evolution, which encompasses the formation (due to 
ambipolar diffusion) and contraction of magnetically supercritical cores in magnetically 
subcritical clouds, we use the models described in CM93, CM94, and CM95. Further details 
are given in Appendix A. Here we summarize key features of these models. 

Clouds are modeled as axisymmetric thin disks, with half-thickness Z(r,t), embedded 
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in a constant-pressure external medium, with the axis of symmetry aligned with the z axis 
of a cylindrical polar coordinate system (r,(f>,z); exact balance between thermal-pressure 
and gravitational forces along magnetic field lines is assumed to hold at all times. These 
simplifying assumptions are based on the results of Fiedler & Mouschovias (1992, 1993), who 
found that initially uniform, spherical or cylindrically symmetric, self-gravitating magnetic 
molecular clouds rapidly flatten and establish force balance (between thermal-pressure and 
gravitational forces) along magnetic field lines. In their typical models, balance of forces 
along field lines was maintained even after the onset of dynamical contraction perpendicular 
to the field lines. The collisional friction of (charged and neutral) grains, which in certain 
cases can be significant (CM94, CM95), is ignored in this paper. We also neglect the effect 
of rotation and magnetic braking (e.g., Basu & Mouschovias 1994, 1995a, b); they will be 
considered in a later paper. 

Model clouds are initially in magnetohydrostatic equilibrium and would remain so 
indefinitely in the absence of ambipolar diffusion. Hence, any evolution is entirely the 
result of ambipolar diffusion, which is "turned on" at time t — 0. The evolution of a 
magnetically supercritical core is followed up to a time t; nt , at which the central density 
exceeds 10 11 cm~ 3 . At this point the pre-PMF calculation is halted. We then proceed as 
described in the next subsection. 



2.2. PMF Phase and Subsequent Evolution; Modification of the Ambipolar 

Diffusion Models 

To calculate this phase of the evolution we use as initial data the values of the physical 
quantities (column density, magnetic field, neutral infall speed, etc.) of the pre-PMF 
epoch at time t- mt . To numerically calculate through the PMF phase of evolution, we use 
the "central sink cell" method originally employed by Boss & Black (1982) to model the 
collapse of nonmagnetic, isothermal clouds. This method was also used in the collapse 
studies of Foster & Chevalier (1993) and Tomisaka (1996). In this method we make the 
central cell of our (stationary) computational mesh a sink cell, with size equal to the mesh 
inner boundary radius r = ri nner : the central sink cell is essentially a "hole" at the origin of 
our computational mesh. Use of the central sink cell keeps the numerical time step St nmn 
[~ r inncr/|' u n( r inner)|, where v n is the infall speed of the neutrals] finite throughout the PMF 
epoch. 

We fix the value of r inner at the radius at which the governing equations of a model cloud 
no longer remain valid. Two fundamental assumptions of our models are isothermality and 
the freezing of magnetic flux into the ions. The isothermal approximation breaks down for 
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n Q ^ 10 11 cm -3 (Gaustad 1963; Hayashi 1966; Larson 1969; Winkler & Neuman 1980); at 
densities above this value, radiative heating of the gas can no longer be neglected. Freezing 
of the magnetic field in the ions is no longer valid when the dimensionless ion magnetic 
attachment parameter r\ = cJiT in , also known as the ion Hall parameter, is ^ 1 (see eqs. [38] 
and [39] of CM93), where uj- x is the ion gyrofrequency and r in is the ion-neutral collisional 
timescale. The CM94 collapse model that neglects the collisional effect of grains (see their 
Figs. 2a and 2f) yields the approximate scaling B zeai ~ 3.5 x 10 3 (n n /2.6 x 10 9 cm~ 3 )°- 4 /xG 
for the equatorial (z = 0) magnetic field. Inserting this relation into equation (40) of CM93 
gives r ; <z 1 for n n ^5x 10 10 cm~ 3 . Thus, co incidentally, the assumption of freezing of 
magnetic flux in the ions breaks down at effectively the same density as the one where the 
isothermality approximation also becomes invalid. We take ri nner to be the radius at which 
n n (r iimor ) ~ 10 11 cm~ 3 . Typically, this corresponds to r inner ^ 5 AU (see § 3.3). 

Mass and magnetic flux accumulate in the central sink from our active computational 
region, consisting of all the cells located at r > ri nner ; the rate at which mass is advected 
into the central sink is ^ 

Mcent = = -2nra n v n \ r=rinnai , (2) 

and the rate at which magnetic flux flows into the central cell is 

$B,ccnt = d^ cent = -2nrB zm v,\ r=n ^ , (3) 

where M cent (t) and $s, C ent(^) are the central mass and magnetic flux, a n (r,t) and B Zfi(l {r,t) 
are the column density and magnetic field strength (in the equatorial plane of the disk), 
and v n (r, t) and t>i(r, t) are the neutral and ion infall speeds. 

Despite the fact that the mass and flux contained within the central sink cell are 
exterior to the active computational region in our models, they still affect the evolution of 
the collapse through their gravitational and magnetic influence on the matter at r > ri nner . 
This has to be incorporated into the calculation of the radial gravitational and magnetic 



3 When I\ < 1 the assumption of flux freezing in the ions is no longer valid. However, one can still consider 
the magnetic flux to be frozen into the electrons, if the electrons are the main charge carriers (as opposed to 
charged grains) and the electron magnetic attachment parameter r c = uj c T cn > 1 (for a discussion, see the 
Appendix of Konigl 1989). At even higher densities than we consider here, T e becomes < 1 and the effect of 
Ohmic dissipation has to be included in the magnetic induction equation (Spitzer 1963; Pncuman & Mitchell 
1965; Norman & Heyvaerts 1985; Nakano & Umebayashi 1986a, b). For detailed discussions of the effect of 
a finite conductivity on the MHD equations for multicomponent plasmas at the high densities encountered 
in collapsing cores and/or protostellar accretion disks we refer the reader to Nakano (1984), Draine (1986), 
Konigl (1989), Wardle & Konigl (1993), Mouschovias (1996), and Wardle & Ng (1998). 
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field components. The radial gravitational acceleration in a thin disk is (see CM93, eqs. 
[31a] and [31b]) 

A roo roo roc roc 

g r (r) = 2iiG— dr' V a n (r') dk J (kr) J (kr f ) = -2nG dk J x {kr) dr' r'a n (r')J (kr') 
dr Jo Jo Jo Jo 

(4) 

where Jo and J\ are Bessel's functions of order and 1, respectively. One can show by 
direct integration of equation (|j) that, for cr n (r) oc r s,T , where s a is a constant between 
-1/2 (corresponding to free-fall during the post-PMF phase; see § 3.3) and -1 (typical of 
the late stages of the pre-PMF phase, see § 3.2), g r (r) = —TaM{r)jr 2 ^ where Mir) is 
the mass enclosed within the radius r and T a is a constant between 0.68 and 1. As an 
example, column density profiles cr n (r) oc r -0 95 typically developed in the inner regions 
of the contracting cores of the models presented by CM94 and CM95 (see also § 3.2), 
corresponding to T a ~ 0.82 in the above expression. Considering the alternative case of a 
highly localized source, such as a uniform-density sphere of radius i?* (<g r inncr ) and mass 
M*, the column density obeys the relation 

n i/2 

, r < , (5a) 



3 M* 



r 

= , r>R*, (5b) 

and one finds from the solution of equation (ffl) that, for r > .R*, g r (r) = —GM+/r 2 , as 
expected. Hence g r {r > r nmcvi t) ~ —GM(r,t)/r 2 for the cases of a power-law column 
density profile and of a localized mass source. As these two extremes tend to bracket the 
range of possible behaviors of cr n (r, t) in the innermost flux tubes of the core, we may use 
the expression 

= _^pt + 2jrG d /- r , 7 7 

r z dr jo Jo 

during the PMF epoch. Extension of the integral in equation (^) to r = is obtained, 
with negligible error, by setting cr n (r < r inn0I .) = a , where o"o is a constant that satisfies 
the condition < a < cr n (r inner ). A description of the numerical method used to solve the 
integral term on the r.h.s. of equation (§) is given in § 2.4 of Morton, Mouschovias, & 
Ciolek (1994, hereafter MMC94). 

As shown in § 3.1.1 of CM93, the r component of the magnetic field at the surface of 
the disk, which appears in the restoring magnetic tension force (see CM93, eq. [28c]), is 
given by 

d z" 00 /"°° 
B r ,z{r) = -~r dr'r'[B zm (r')-B Ie{ ] dk J (kr)J (kr') 
dr Jo Jo 

= / dkJ^kr) dr'r'[B zm {r')-B Tei ] J (kr') , (7) 
jo jo 
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where the constant _B re f is the external (background) magnetic field strength at r, z — > oo. 
Typically, B ZCOL {r) ^> B Tcl in the inner flux tubes of a contracting core (CM94, CM95; Basu 
& Mouschovias 1994, 1995a, b). Q 

By analogy with the discussion of the gravitational field g r (r) in the preceding 
paragraph, calculation of B r ^z{r) for the bracketing source magnetic field profiles 



B zm (r) oc r' 
for r < i?*, and 5^ ieq (r 



Sb , and that of a central point-flux source [i.e., B 



B+ 



const 



■ for r > i?*], equation (|7D yields, for both cases, 
B r ^z{r) ~ $s(r)/27rr 2 , where $s(r) is the total magnetic flux enclosed within radius r. It 
follows then, using a derivation similar to that of equation that B r Z {r > r inncr ,t) after 
PMF can be written as 



BrAr) 



_B,cent 



2nr< 



— / dr' r' [B zm (r') - B rc{ ] / dk J (kr) J (kr') . (8) 
dr Jo Jo 



Extension of the integral on the r.h.s. of equation (H) to r < dinner is done in a fashion 
similar to that discussed above in connection with equation (f|). 

The effect of the gravitational field of the central point mass also needs to be included 
in the equation of quasistatic equilibrium along magnetic flux tubes. Assuming balance of 
gravitational and thermal-pressure forces in the z direction within the disk yields 



-Pext + 



Biz 



•57T 




" ri 2 , GrM cen tp n 

+ -Gff n H 



7T 
— I 

2 



<7 r 
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(9) 



In deriving equation (^|) we have used, as in CM93-CM95, the "one-zone 
approximation" for the relation between the mass density p n and the column density 
cx n : a n (r,t) = f^zfy p n ( r ,t)dz = 2p Q (r,t)Z(r,t). The first three terms on the r.h.s. of 
equation @, present even in the absence of a central point mass (see CM93, eq. [26]), are, 
respectively, the constant external pressure, the magnetic squeezing associated with the 
radial field component at the disk surface, and the self-gravitational stress of the matter 
contained in a flux tube; the last term on the r.h.s. of equation @ is the tidal gravitational 
stress corresponding to the z component of the gravitational field of the central point mass. 
Equation (0) simply states that the pressure in the equatorial plane is equal to the sum of 



4 Note that, under the conditions B 2jeq (r) 3> B le { and B z ^ eq (r) oc cr n (r), which in our model calculations 
are satisfied before PMF for r r COTC , equations (||) and (0) yield B ri z(r) oc g r (r). The condition 
-Bz,eq(f) c n (r) expresses flux freezing into the neutrals, and the relation B r z oc g r was derived under 
this assumption by Li & Shu (1997) and Zweibel & Lovelace (1997), who considered model clouds with 
B IC { = 0. This relation was also employed by Basu (1997) in modeling pre-PMF core collapse. 



-12- 



the external thermal plus magnetic stresses and the total gravitational stress acting on each 
flux tube; this equation is used to calculate p n (r,t). f] 

As it turns out, the first two terms on the r.h.s. of equation @ have a negligible 
effect on the disk evolution in our simulations. In particular, the ratio of P ex t and the 
self-gravitational stress is generally <C 1 in the inner flux tubes of the core (see eq. [A2]). 
Furthermore, we have verified that, with the exception of the innermost zones of the 
active computational grid near the end of our typical simulation, the magnetic squeezing 
term is smaller than the total gravitational stress in equation (H). Since the properties 
of these zones have only a small effect on the core structure at that epoch (in particular, 
the ambipolar diffusion-driven shock that propagates through the core has by that time 
reached much larger radii; see §3.3), we have, for simplicity, omitted this term altogether 
in our calculations. As we show in Appendix B, when the core approaches free-fall collapse 
under these conditions, equation (||) implies p n oc r~ 2 . This behavior, which differs from the 
dependence p n oc r~ 3 / 2 that characterizes spherical infall onto a point mass, results from 
our use of the one-zone approximation in the relation between the column density and mass 
density for the thin-disk cloud model. 

A final modification that is required for the post-PMF phase of collapse is in the ion 
force equation, from which the ion-neutral drift speed v u = v\ — v n is derived (see CM93, 
eqs. [50] and [51]). This has to do with the fact that, as first pointed out by Mouschovias 
& Paleologou (1981), the ion-neutral collision rate (aw) n { = 1.7 x 10 _9 cm 3 s _1 = const 
(used in the ambipolar diffusion models of CM93-CM95), which was derived by using the 
well-known Langevin approximation (see, e.g., Gioumousis & Stevenson 1958; McDaniel & 
Mason 1973), is valid only so long as the drift speed satisfies the criterion (<jw) ni /\v D \ > cr geo , 
where <7 geo is the geometric cross section for ion-H 2 collisions (= 7r[a ion + i, where 
a ion and a Ha are the ion and H 2 -molecule radii, respectively). Therefore, for drift speeds 
v D > v Dctit = (cw) n i/cr gco , the quantity (7g e o|f D | must be used for the collision rate, and 



5 For computational convenience, we continue to use the quasistatic approximation for motions in the z 
direction even when there is dynamical collapse in the r direction. To test the validity of this assumption, we 
have also run models that did not employ the quasistatic assumption along flux tubes and instead allowed 
for acceleration of the neutrals in the z direction. In these models there was initially a short-lived phase of 
vertical dynamical relaxation and oscillation in the central flux tubes, followed by rapid reestablishment of 
balance of forces along flux tubes, in agreement with equation (|^). The overall evolution of these models 
differed only slightly from those in which balance of vertical forces was assumed at all times. However, the 
computing time needed to run models that did not assume vertical force equilibrium was typically an order 
of magnitude longer than for models that always used the quasistatic approximation. 
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instead of using equation (51) of CM93, the appropriate equation for the drift speed is 



1.4 



(m„ 12 a.m.u.) 

1 + 0.067^-^ r 

(mi/30 a.m.u.) 



mag,r 



1/2 



t)a n (r,t) 



(10) 



where F mag r (r,t) is the total (radial) magnetic force per unit area, which includes both the 
magnetic tension and pressure stresses (see eq. [28c] of CM93). The factor 1.4 in equation 
(|I0|) reflects the fact that we are neglecting the inertial effect of a 20% He abundance in the 
neutral-ion collision timescale r n i (see eq. |TJ). For the various species of ions we include 
in our models (see Appendix A) we adopt a "generic" ion radius of 0.8 A; this yields a 
critical drift speed f Dcrit = 12 km s _1 , only slightly different from the value of 10 km s _1 
obtained by Mouschovias & Paleologou (1981), who had assumed a plasma consisting solely 
of Na + ions. Hence, when v D (r, t) > f Dcrit , we use equation fllCf) ; otherwise, equation (51) of 
CM93 is used. As we show in § 3.3, ambipolar diffusion can continue to be effective even in 
dynamically contracting cores, and drift speeds v D > v D crit may be attained in certain core 
regions during the post-PMF epoch. 

Spatial discretization and time integration of the equations governing the evolution of 
a model cloud are carried out by the methods described in MMC94 — with the exception 
that spatial derivatives in the first computational cell, with inner boundary r = r inner , are 
calculated by using one-sided differences instead of the three-point technique described in 
§ A2.4 of Mouschovias & Morton (1991). All computations were performed on an SGI 
R-4000 Indigo workstation; running in background, a typical model calculation took ~ 2 
weeks to form a central protostar with a mass of 1 M & . 

As we show in § 3.3, ambipolar diffusion significantly alters the evolution of a collapsing 
core following PMF. To verify that the results of our typical simulation represent a real 
effect that arises from a physical diffusion process and are not an artifact of numerical 
diffusivity, we have run a parallel model that corresponds to flux freezing into the neutrals 
for all times At > 0. (This was accomplished simply by setting v- x = v Q in our numerical 
code.) The results of this calculation (detailed in § 3.3) have revealed no evidence of 
numerical diffusion arising from our finite discretization scheme: each neutral fluid element 
was found to evolve with constant mass and magnetic flux throughout the entire simulation, 
as would be expected for a fluid with a frozen-in magnetic field. In § 3.3 we therefore use 
the results of the frozen-flux cloud calculation as an aid in isolating and interpreting the 
effects of ambipolar diffusion in our typical model (see also Appendix C). We note in this 
connection that we have also run other ambipolar-diffusion models with different mesh sizes 
and core resolutions and obtained essentially the same results as in our typical model. This 
further confirms the absence of numerical-diffusion effects in our simulations. In addition, 
we have tested the effect of varying the assumed profile of B z eq for r < r inner in calculating 
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B r z from equation (S), and confirmed that our quantitative results were insensitive to the 
detailed form of that profile. 

3. Numerical Simulation 

3.1. Model Parameters 

We present a typical model of the dynamical evolution of a magnetically supercritical 
core through the PMF phase. The values of the dimensionless parameters of this model are 
listed in Table 1 and Table 2. They are similar to that of model Auv presented in CM95. 
The values of the free parameters selected could be used to represent a molecular cloud 
with temperature T = 10 K, radius Rq = 4.29 pc, and total mass M d = 98.3M Q . For a gas 
of H2 with a 20% He abundance the mean mass of a neutral particle is m n = 2.33 a.m.u. 
and the isothermal speed of sound is C = (/c B T/m n ) 1//2 = 0.19 km s _1 . The initial central 
density n niC0 , column density o" niC0 , and magnetic field strength -B 2ieq ,co are 2.6 x 10 3 cm -3 , 
5.59 x 10~ 3 g cm -2 , and 35.3 fiG, respectively. The dimensional values of the parameters 
used in the calculation of the equilibrium abundances of charged particles, such as the 
abundances of the different atomic species, cosmic-ray and UV ionization rates, chemical 
and charge-transfer reaction rates, etc., are the same as those cited in § 3.1 of CM95, with 
the exception that the probability V\ of ions sticking onto grains (see CM93, eqs. [57b] and 
[57c]) has been changed from 0.9 to 0.99. For the purposes of calculating ion abundances, 
we assume a uniform population of grains with radius a = 3.75 x 10~ 6 cm. [] 

3.2. Pre-PMF Phase: Supercritical Core Formation and Contraction 

Models of the self-initiated (due to ambipolar diffusion) formation and contraction of 
magnetically supercritical cores have been presented and discussed at length by Fiedler & 
Mouschovias (1993), CM94 and CM95, and Basu & Mouschovias (1994, 1995a, b); we refer 
the reader to these papers for more detailed descriptions of the physics of the formation 
of protostellar cores in magnetically supported interstellar clouds. During the pre-PMF 
phase, the evolution of the core of the typical model cloud presented here is similar to that 



6 Although we do not account for grain-neutral friction in this model, our results would be little changed 
if we had instead used grains with radius a J> 7 x 10~ 6 cm and included the effect of grain-neutral collisions, 
as grains of this size or larger contribute only marginally to the total collisional force on the neutrals (see 
model 2 of Ciolek 1993). 
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of model A uv discussed at length in § 3.2.1 of CM95; for brevity, we summarize here only 
some of the features of the core in the pre-PMF phase. 

Figures la-li display, respectively, the density (normalized to ra n ,co), column density 
(normalized to <7 niC o), mass-to-flux ratio (in units of the critical value for collapse), magnetic 
field strength in the equatorial plane of the disk (normalized to B zcqjC0 ), r component of 
the magnetic field at the surface of the disk (normalized to -B z , e q,co)> infall speed of the 
neutrals (normalized to C), drift speed (normalized to C), mass infall rate M (in units 
of M© Myr~ x = 0.65C 3 /G), and ratio of local cloud vertical half-thickness Z and radius 
r as functions of r/Ro, at eleven different times tj. The times tj correspond to when the 
central density n njC (tj) = 10% niC0 ; dimensionally, these times are t — 0, ti — 7.56 Myr, 
t 2 = 9.29 Myr, t 3 = 9.56 Myr, t 4 = 9.607 Myr, t 5 = 9.6167 Myr, t 6 = 9.6169 Myr, 
t 7 = 9.6190 Myr, t 8 = 9.61977 Myr, tg = 9.61981 Myr, and t w = 9.61983 Myr. An asterisk 
on a curve in Figure 1 locates the instantaneous radius R C rit{tj) inside which the total- mass 
to total-flux ratio is equal to the critical value for collapse, i.e., 

(M\ 1 



- 



An open circle locates the instantaneous critical thermal (~ Jeans) lengthscale AT, C rit(^) 
[= C 2 /2Go~ UtC (tj), where cr n ,c(tj) is the central column density at time tj] Mouschovias 1991]. 

During the later stages of the evolution (t ^ ti), the core develops a uniform central 
region [of radius ~ AT,crit which is continually shrinking in both size and mass, and 
a "tail" of matter and magnetic field left behind by the central region. Inside the "tail" 
region, near power-law behavior emerges, with rflnp n /<ilnr —2, d\na n /d\nr pa —1, 
and d In B z>cq /d In r pa —1 deep inside the core, and with d\np n /d\nr pa —1.5, 
d\na n /d\nr pa —0.7, and d In B zeq /d In r ~ —0.6 further out, where most of the core mass 
is contained. As noted by Basu (1997), for t ^ t±, the column density and z component of 
the magnetic field within the inner part of the core (r i? CT it) are well approximated by 
the relations 



Q~n,c(t) 

(l + [r/2A T , cr (t)] 2 ) 

and 

B z ,eq,c \t) 

(l + [r/2A T)C r(0] 

where a njC (t) and B zeqc (t) are the values of a n and B z eQi in the uniform central region 
at time t (see Figs. 16 and Id). PMF occurs when At, ct (^) ~~ * 0. In the inner core the 
magnitude of the retarding magnetic force decreases relative to the gravitational force, 



*n(r,0-- nxl/2 (12) 



U t) * - "''^^ , (13) 
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becoming comparable to that of the thermal-pressure force. Further out, at radii ~ Rcnt(tj), 
the magnetic force is the primary source of support against self-gravity. 

By the time t\o the core has a radius -R cr it(^io) = 2.59 x lCT 2 i?o = 0.111 pc and a mass 
M core = 5.37 M Q . These values are in excellent agreement with typical observations of dense 
ammonia cores in star-forming molecular clouds (e.g., Ladd, Myers, & Goodman 1994). 
Separation of the magnetically supercritical core from the massive, subcritical envelope is 
demonstrated in Figure lh, which shows that the mass infall rate M is reduced by more 
than four orders of magnitude for r > R CT it(tj). This is due to the effective freezing of the 
magnetic field in the neutrals at large radii, which is a consequence of the comparatively 
high degree of ionization brought about by the penetration of the external UV radiation 
field into the optically thin cloud envelope (see CM95 for a discussion). [] 

We note the following results of interest for the pre-PMF phase of evolution. First, as 
can be seen from Figures 1/ and lg, the magnitude of the infall speed of the neutrals in the 
inner regions of the collapsing core is comparable to the isothermal speed of sound C (and 
thus also to the fast-magnetosonic speed v ms = [C 2 + v 2 p] 1 / 2 , where va = -B^eq/^vrpn] 1 / 2 is 
the Alfven speed, since v ms is nearly equal to C in these regions). Moreover, at time ti , the 
inward acceleration a n of the neutrals lies in the range 0.25g r — 0.5g r in the inner regions 
of the core, and \v n \ is a fraction ~ 0.55 of the free-fall speed (~ 1.98 C) at r = ri nner (see 
§3.3). Hence the core is indeed dynamically collapsing (although not freely falling), so, 
as we noted in § 1, treating a supercritical core as being in near hydrostatic equilibrium 
at PMF is not a valid approximation for molecular cloud cores that form and evolve by 
ambipolar diffusion. 

Second, we note the increase of the drift speed v D in the inner flux tubes for t ^>t 6 (see 
Fig. lh), which indicates that ambipolar diffusion can continue to operate even during the 
dynamical collapse of the core. This is primarily due to the phenomenon of ion depletion, in 
which ions become increasingly attached onto grains at higher densities, thereby reducing 
the relative abundance of ions in the gas phase. This, in turn, results in a higher diffusion 
rate (see CM93, CM94) and leads to an increase (for i > t 6 ) of the mass-to-flux ratio in the 
flux tubes that thread the inner core (see Fig. lc). (The effect of ambipolar diffusion on 
the mass-to-flux ratio during the approach to PMF is also discussed in Basu 1997.) As we 
show in the next subsection (see also § 1), ambipolar diffusion in the interior flux tubes is 



7 The outer envelope may be supported even better than we have calculated in view of the fact that the 
UV photoionization rates could be higher than those listed in Table 2, in which case the fractional ionization 
values in the outer regions of the cloud (corresponding to visual extinctions <^ 4) would be slightly larger 
(Ruffle et al. 1998). 
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increased dramatically after PMF and plays a significant dynamical role in the subsequent 
evolution of the core. 

Finally, we notice that by the time t\o, Z/r < 1 throughout the whole cloud, except 
for a very small inner region of the core (r/Ro <^ 10~ 7 ; see Fig. li). Note that Z/r > 1 
does not necessarily indicate an inconsistency in our use of the thin-disk approximation 
to model clouds. As discussed in CM93 and CM95, this approximation will be valid 
so long as any scalar quantity f(r) (such as n a , a n , B z eq , etc.) satisfies the condition 
f(r,t)/\df/dr\ > Z(r,t). Examination of Figures la, 16, and Id reveals that this condition 
is satisfied even in the uniform-density central region because df/dr m for r — > before 
PMF. As it turns out, the innermost region of the core, which at later times happens to 
contain densities n n > 10 11 cm -3 (see Fig. la), is in any case excluded from our calculation 
of the post-PMF core evolution for the reasons outlined in § 2.2. 

3.3. Point-Mass Formation and Protostellar Accretion Phase of Core 

Evolution 

The physical data corresponding to the model at time t int = t w are used as input for 
the protostellar-accretion phase of the calculation. (Our results for the post-PMF epoch are 
relatively insensitive to the particular value of tint we choose; for instance, we find that the 
resulting evolution differs only marginally if we use t- m t = ts instead.) The boundary of the 
central sink cell is taken to be r inner = 8.24 x 10 _6 i?o = 7.31 AU. (The post-PMF evolution 
is also relatively insensitive to the value of r inner so long as r inner is much smaller than the 
core radius R crit ~ 0.1 pc; for instance, the results for r inner ~ 20 AU are very similar to 
those of our typical run.) The density at this position and time is 5.0 x 10 10 cm -3 (see Fig. 
la), whereas the mass contained in the central cell at this time is M cen t(^io) = 8-5 X 10~ 4 M Q , 
which is negligible in comparison with the total core mass M core = 5.37 M Q . Similarly, 
^B,cent(^io) = 6.72 x 10 26 Mx, which is several orders of magnitude smaller than the total 
core magnetic flux $s iCO re = 1-72 x 10 31 Mx. 

The evolution of the mass in the central sink cell is shown in Figure 2a as a function of 
At = t — t w . Point-mass formation, signaled by rapid growth of M cent with time, takes place 
at At < 200 yr. R - - because of the nonzero value of r inner , we can only place an upper 
limit on this time. We terminate the simulation at At = 1.53 x 10 5 yr, when M cent = 1M . 



8 Therefore, for At ;> 10 3 yr, the quantity At can be considered to be the time elapsed since PMF. This 
makes our time unit At the same as the time t typically used in previously published self-similar collapse 
models (see references cited in § 1), which set t = at PMF. 
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The accretion rate M cent onto the protostar (in units of M Q Myr -1 ) is shown in Figure 
26 as a function of At (solid curve). As in the previous (nonmagnetic and magnetic) 
simulations of dynamical collapse cited in § 1, M cent rapidly increases imediately after PMF 
due to the enhancement of the gravitational acceleration brought about by the appearance 
of the protostar. M cent reaches a maximum value ~ 9.4 M Q Myr -1 , which is less than the 
theoretical upper limit of 13C 3 /G = 20 M Myr -1 estimated by Basu (1997) for collapsing 
cores in which ambipolar diffusion remains operative. M cent decreases for At ^ 4 x 10 3 yr. 
As we show below, the accretion rate during this time is strongly affected by an ambipolar 
diffusion-induced hydromagnetic shock that propagates outward and slows the infall of 
matter at larger radii. We also plot in Figure 26 the central accretion rate of a collapse 
model (dashed curve) that began with the same initial data as our typical model, but with 
v D set equal to zero for At > (i.e., for this model, as discussed in § 2.2, the magnetic 
flux was artificially forced to remain frozen into the neutrals during the post-PMF epoch). 
Finally, in Figure 2c, we plot M ccnt (in the same units as in Fig. 26) as a function of 
M cent /M . The dashed curve again displays the accretion rate for the frozen-flux model. 

Comparing the two models shown in Figures 26 and 2c, we note that M cent is the same 
for At < 10 3 yr: this has to do with the fact that the ambipolar diffusion timescale r AD 
in the typical model is much longer than the gravitational contraction timescale r gr during 
this period (see Figs. 3a and 36), so that diffusion does not greatly alter the evolution. 
However, for At > 10 3 yr, r AD ~ r gr , so ambipolar diffusion is much more effective and leads 
to a redistribution of the mass in the inner flux tubes of the core (see below). This, in turn, 
results in a reduction of the mass accretion rate onto the central protostar. Specifically, 
by the end of the simulation, M cent in the ambipolar-diffusion model is ~ 60% of the 
accretion rate in the frozen-flux case. The ambipolar diffusion-induced decrease in M ccnt 
comes on top of the monotonic decline exhibited in the flux-frozen case for At ^ 10 4 yr. 
The behavior of the frozen-flux model is consistent with the results of other nonmagnetic 
and magnetic collapse calculations (e.g., Hunter 1977; Foster & Chevalier 1993; Tomisaka 
1996; Safier et al. 1997; Li 1998). As discussed by Basu (1997, § 5), the post-PMF decline 
in M ccnt is due to the fact that, at the time of PMF, the outer mass shells of the core are 
not as strongly accelerated inward (and are therefore moving more slowly than the inner 
mass shells) because of their much larger initial distance from the central point mass. The 
maximum values of M cent in both of these numerical models agree with those found in the 
self-similar magnetic collapse solutions presented in Contopoulos, Ciolek, & Konigl (1997, 
hereafter referred to as CCK97). 

We now show that, during the post-PMF epoch, ambipolar diffusion operates effectively 
in the inner core even as it continues to undergo dynamical collapse. As we did in §1, 
we again examine the ambipolar diffusion timescale r AD = r/v D , which we now express in 
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a manner similar to that of Mouschovias (1989, 1991). The force (per unit area) on the 
neutrals is 

2 da n a n 

a n a n = a n g r - C — H v D , (14) 

OT T n [ 

where a n = dv n /dt is the acceleration of the neutrals. The first term on the r.h.s. of 
equation ( jTjj ) is the gravitational force (including both self-gravity of the gas and the 
gravitational field of a central point mass), the second term is the thermal-pressure force, 
and the last term is the frictional force due to collisions between neutrals and ions. Dividing 
this equation by r, one obtains 



(T gr /r n i) 



T'gr/'ace) (''"gr/'s) 



(15) 



where r gr (= [r / \g r \\ 1 ^ 2 ) is the gravitational contraction timescale (referred to as the 
dynamical timescale in CM94 and CM95), r acc (= [r/ |a n |] 1//2 ) is the acceleration timescale, 
and r s (= \{C 2 /r) \dhi<j n / dr\]~ l l 2 ) is the sound crossing time. Equation ( |15"D is the same 
as equations (19) and (20) of Mouschovias (1989), except that we use the gravitational 
contraction timescale r gr instead of the free-fall time r ff . Now, during the later stages of 
dynamical collapse, the ratio r gr /r acc in the denominator approaches unity and tends to 
increase T AD /r gr (the ratio [T gr /r s ] 2 is negligible). However, effective diffusion (corresponding 
t° T Ao/ T gr ^ 1) can s ^ m occur if the ratio r gr /r ni in the numerator of equation (|T5| ) is 
small during this period — due either to a decrease in r gr or to an increase in r n j. Hence, 
if r gr /r n i <C 1, ambipolar diffusion can still redistribute mass and flux in a dynamically 
collapsing core. | That this is indeed the case in our collapse model can be seen in Figures 
3a and 36, which show, as functions of At, T AD /r gr , r gr /r ni , and the total (nondimensional) 
magnetic flux contained within cells 1 and 5 of our computational mesh. Examination of 
these figures reveals that r gr /r ni decreases with increasing At. This is caused solely by the 
decrease in the value of r gr oc |g r |~ 1//2 oc M~H 2 (see eq. 0). On the other hand, r ni oc n^ 1 
(see eq. |TJ) changes very little during this period due to the fact that, as noted in § 3.2, 
ni ~ const for n n ^> 10 7 cm~ 3 . The rapid increase of M cent with At following PMF (see Fig. 
2a) thus leads to a strong decrease in r gr /r ni , resulting in a "revitalization" of ambipolar 
diffusion in the inner flux tubes during this phase. For At ^ 10 3 yr, T AD /r gr ~ 1, and the 
growth of the magnetic flux contained within these two cells is halted. The flux then begins 
to "pile up" in the inner cells of the computational mesh. (The revitalization of ambipolar 



9 The criterion r gr /r n i <C 1 is equivalent to requiring Am,ct/^a -C 1, where Am,ct ~ vat {{ ~ fA^gr is the 
critical magnetic lengthscale, and Aa ~ «AT n i is the AlfVen lengthscale, as defined by Mouschovias (1991). 
The condition Am,ci- *C Aa for effective ambipolar diffusion during rapid core collapse was put forth by 
Mouschovias (ib., § 4.1). 
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diffusion after PMF can also be understood in terms of the scaling of the magnetic field 
following point mass formation; see Appendix C.) 

The effect that the piled-up magnetic flux has on the inner flux tubes of the typical 
model can be seen in Figures 4a and 46, which display, respectively, B Zfi J B z ^ c0 and 
On/on,co as functions of At for cells 2, 5, 10, 15, 20, 25, and 27 of the computational mesh. 
Figure 4a shows that the piling up of magnetic flux begins close to the axis of symmetry 
and then, with increasing time, becomes noticeable at progressively larger radii. Hence, 
after PMF, ambipolar diffusion causes a front of magnetic flux to move outward from 
the innermost flux tubes. Ambipolar diffusion is so effective behind this hydromagnetic 
disturbance (HMD) that the inward advection of magnetic flux stops — with the field 
lines being held in place or else moving slowly outward, in either case with \v{\ <C |i> n |. 
This behavior is reproduced by the self-similar collapse solutions of CCK97, which also 
incorporate the effect of ambipolar diffusion on the post-PMF core evolution. The magnetic 
field strength B ZjSq is greatly enhanced after the passage of the HMD in each cell. The 
effect of the HMD on the behavior of the column density varies with the distance from the 
center. In particular, the column density in cells 2 and 5 shows no effect from the passage 
of the magnetic field front (see Figs. 4a and 46). This is because ambipolar diffusion is 
so rapid in these cells that very little of the increased magnetic force from the enhanced 
field is transmitted to the neutrals through ion-neutral collisions. However, after the HMD 
reaches larger distances, the collisional coupling between the two fluids increases (see Fig. 
Qh). Furthermore, the piling- up of flux increases the local field strength at these distances 
by a larger relative factor. As a result of these two effects, the increased magnetic force 
acting on the neutrals is able to temporarily decelerate the infalling matter (from slightly 
supersonic to subsonic infall speeds; see Fig. 6e) for cells with index I > 8. For these cells, 
er n increases immediately after the passage of the front. The effect of the HMD at larger 
radii can be described in terms of a hydromagnetic shock that propagates in the weakly 
ionized gas. (As noted in § 1, the formation of such a shock was predicted by Li & McKee 
1996; we further analyze the nature of this C-type shock in § 4.2.) After traversing the 
shock, the neutrals are gradually reaccelerated and eventually reach free-fall speeds. /,From 
Figure 46 we find that s At = dlncr n /dln At lies in the range —0.76 ^ s At ^ —0.55, which 
can be compared with the value s At = —1/2 that characterizes self-similar collapse models 
of thin, axisymmetric disks (Li k Shu 1997; CCK97). 

The position of the hydromagnetic disturbance r HMD (t) in units of the cloud radius R 
is presented in Figure 5a as a function of At. The speed v HMD (= dr HMD /dt, in units of C) 
of the front in the reference frame of the central protostar is shown in Figure 56. (The time 
intervals between successive data points in these figures are larger than those in Figures 
2-4. Figures 5a and 56 therefore show a time-averaged motion of the HMD and do not 
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exhibit the steady, shock-induced oscillations seen in Figures 2-4 at later times.) Early on, 
the radially outward motion of the front is nonsteady, and depends on the rate at which the 
flux external to the front enters and piles up in the region of rapid ambipolar diffusion. At 
later times, v HMD < C. The Mach number of the hydromagnetic disturbance in the reference 
frame of the upstream gas, .Mhmd = (f HMD — v UtU )/C (where v QtU is the infall speed of the 
neutrals just ahead of the front), is shown in Figure 5c as a function of At. The motion of 
the HMD relative to the upstream gas is supersonic, and, for At ^ 10 4 yr, Mhmd — 2.8. It 
is also of interest to note that, in the vicinity of the front, the fast-magnetosonic speed t> ms 
is typically ^ 2C. Hence the hydromagnetic Mach number -M mSj HMD = (f HMD — v n )/v ms of 
the HMD relative to the upstream neutrals is also > 1. 

Spatial profiles of various quantities in the typical model are shown in Figures 6a-6i 
as functions of r/R , at seven different times At. As noted earlier, At is essentially 
the time elapsed since PMF. For this typical model, At = 0, Ati = 2.17 x 10 2 yr, 
At 2 = 6.11 x 10 2 yr, At 3 = 1.54 x 10 3 yr, At 4 = 3.81 x 10 3 yr, At 5 = 2.38 x 10 4 yr, and 
At 6 = 1.48 x 10 5 yr. An asterisk on a curve locates, as in Figures la-li, the instantaneous 
position of the critical flux tube R CT i t (Atj). The core radius does not change from its 
initial value of R crit {At = 0) = R crit {t w ) = 2.59 x Kr 2 i? = 0.111 pc = 2.30 x 10 4 AU, 
providing further evidence for effective core-envelope separation, as discussed in earlier 
models presented by Mouschovias and coworkers (see references in § 1). 

Figure 6a shows the total mass M in M . For Atj > the curves flatten for small 
radii, revealing point-mass formation since M(r)/r ^ for r — > 0. For r > R cr i t (Atj), M(r) 
changes very little as the magnetically subcritical envelope continues to be supported by 
magnetic forces. There is a small transfer of mass from the envelope to the supercritical 
core during the post-PMF epoch; by the time At 6 the core mass has increased by 8% to 
5.80 M Q . Figure 6b displays B zeq normalized to -B 2iCqiC o, Figure 6c shows the ratio o" n /<T nc0 , 
and Figure 6d depicts B rt z, the r component of the magnetic field at the upper surface of 
the disk, normalized to -B 2 ,eq,co- The infall speed v Q of the neutrals (normalized to C) is 
shown in Figure 6e, the ion-neutral drift speed v D (also normalized to C) is displayed in 
Figure 6/, and the mass accretion rate M in M Myr^ 1 (= 0.65C 3 /G for the temperature 
and gas composition assumed in this model) is exhibited in Figure 6g. The ratio r n i/r gr is 
plotted in Figure 6h, and the ratio Z/r is shown in Figure 6i. 

Taken together, Figures 6b-6i show the effect of ambipolar diffusion and the outward 
progression of the HMD (and the resulting hydromagnetic shock) on the evolution of 
the core. For At < At 2 , ambipolar diffusion has not yet been revitalized, and mass and 
magnetic flux are continually advected inward. During this period, a n and B z ^ eq in the 
core decrease with increasing At. B r z increases as field lines are bent inward while being 
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carried along with the neutrals. For Atj > At 2 , ambipolar diffusion has become rapid 
enough to halt the advection in the innermost cells. Subsequently the front of piled-up 
magnetic flux expands outward to larger radii, dramatically enhancing the value of B zeq 
in its wake. Bending of the interior field lines is slowed and temporarily arrested in the 
region of rapid ambipolar diffusion behind the HMD for At 2 <^ At <; At*, (the behavior for 
At > At 5 is considered below). Further out from the center, the magnitude of the local 
gravitational field strength decreases while the degree of ionization increases, resulting in 
less efficient ambipolar diffusion and more effective collisional coupling between the ion and 
neutral fluids. At times At ^ At 4 the HMD is propagating out into the regions of the core 
( r /Ro ^ 1CT 4 ) where r ni /r gr ^ 1 (see Fig. 6h). At these radii both the collisional coupling 
of the ion and neutral fluids and the enhanced magnetic field strength behind the HMD 
are large enough to affect the inflow of the neutrals, and a shock forms. Our numerical 
calculation thus confirms the evolutionary sequence originally outlined in Li & McKee 
(1996). A local maximum in cr n occurs because the neutrals are "hung up" (i.e., decelerated; 
see Fig. 6e) in the region of stronger magnetic field strength, and their infall is interrupted 
as they are forced to diffuse through the stationary, or slowly expanding, ions and field lines 
contained in the HMD (v D |t> n | and \vi\ <C |i> n ] behind the HMD; see Fig. 6/). After 
the neutrals diffuse through the shock, they resume dynamical collapse toward the origin: 
at time At 6 , |a n | ^ 0.9(% r | for r 4 x 10~ 5 i? = 35.5 AU, a n oc r~ 1/2 , and v n oc r~ 1/2 
(which is the characteristic behavior of free-fall collapse induced by a central point mass, 
also found in self-similar models of gravitational collapse in isothermal disks; Li & Shu 
1997; CCK97). At that time |u n | is a fraction ~ 0.94 of the free-fall speed (~ 76.46 C) at 
f = dinner- Depletion of matter, caused by the interruption in the infall of the neutrals at 
the shock front, is the reason for the local minimum in a n that occurs between the shock 
front and the free-fall region. This causes the accretion rate behind the shock to decrease 
for At ^ At 4 . Hence ambipolar diffusion, through the action of the hydro magnetic shock, 
has increased the magnetic field strength in the collapsing core for r ^ 10~ 4 i? (— 88.7 AU) 
and reduced the accretion rate onto the central protostar — exactly the opposite of what 
one would have expected a priori. 

For comparison we also plot in Figures 6b, 6c, and 6d (dashed curves), for the model 
in which flux freezing into the neutrals was assumed to hold throughout the post-PMF 
epoch, the profiles of o~ n , B Zjeq , and -B r ,z, respectively, at the time At^ the mass of the 
central protostar reached M cent = 1M . (For this model, Ate — 1-23 x 10 5 yr.) In the 
absence of ambipolar diffusion, the magnetic flux in the core is continually advected into 
the central sink, and B Zfiq is greatly reduced from the value that characterizes the core of 
the typical ambipolar-diffusion model (see Fig. 66); B r< z, on the other hand, is much larger 
in the core of the frozen-flux model (see Fig. 6c) because the field lines continue to be 
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bent inward as mass and flux are advected into the central sink. Examination of Figure 
6c reveals that the column density of the ambipolar-diffusion model is significantly larger 
than that of the frozen-flux model for r ^ 1.6 x 10~ 4 i?o — 142 AU, i.e., the region just 
behind the shock front. For r <^ 142 AU the opposite occurs: cr n in the model that includes 
ambipolar diffusion is much smaller than in the frozen- flux case. This is because more of 
the core mass remains further away from the center in the ambipolar-diffusion model, being 
trapped at larger radii by the pressure of the enhanced magnetic field in the wake of the 
HMD. Hence, ambipolar diffusion, by giving rise to a hydromagnetic shock, has significantly 
altered the evolution of the column density profile in the core. Moreover, a hydromagnetic 
shock does not occur in the core of the frozen-flux model. This is contrary to the results 
of the self-similar model of Li & Shu (1997). This discrepancy arises because the initial 
conditions assumed by Li & Shu correspond to a stationary configuration rather than to the 
dynamical one indicated by the numerical simulations of the pre-PMF phase of collapse (see 
§ 3.2). In fact, as discussed in CCK97, shocks also do not form in the absence of ambipolar 
diffusion in self-similar models that do not use the initial singular-isothermal-disk setup 
of Li & Shu. Finally, we note that, for r ^ i? crit (At), the magnetic field is the same in 
both models, reflecting the fact that ambipolar diffusion is so ineffective in the massive, 
subcritical envelope that the field there can be considered to be frozen into the neutrals at 
all times. 

Figure 6b shows that, by the time At 6 , B zeq has been significantly reduced for 
r ^ 2.5 x 10~ 4 -R — 222 AU. This is due to the fact that the drift speed v D (r) exceeded the 
critical value f Dcrit ~ 12 km s _1 in this region of the core for Ats < At < At$. Therefore, 
during this period, the neutral-ion collisional force scaled quadratically with v D rather 



than linearly (compare eqs. [T(| and [0), so the frictional force between the ions and the 
neutrals was significantly increased in this part of the core. This had the effect of refreezing 
the magnetic flux into the neutrals in the innermost computational cells, resulting in the 
field being dragged by the collapsing neutrals into the central sink and therefore decreasing 
in strength in these cells. However, further out (r > 2.5 x 10~ 4 -Ro), where v D (r) < v Dciit , 
the bulk of the flux front remains largely unaffected by the collapse of its "floor." By 
Ai 6 , though, v D has fallen below v Dciit (see Fig. 6/) even in the innermost cells of the 
computational mesh. Bending of the field lines resumes in this region of the core at that 
time, with a subsequent increase in B Ty z (see Fig. 6d). 

We note that beyond the core radius [i.e., for r > R clit (At) ~ 2.52 x 1CT 2 R = 0.11 pc] 
\v n \ <C C (see Fig. 6e), reflecting effective support of the massive, subcritical envelope by 
magnetic forces. Ambipolar diffusion in the envelope is so ineffective (resulting from the 
much greater relative abundance of ions, due to ionization by the external UV radiation 
field) that t) D <C (see Fig. 6/), and the magnetic field there can be considered to be 



-24- 



frozen into the neutrals; furthermore, B r z <C B Zjeq in the envelope (see Figs. 6b and 
6d), indicating that the field lines there remain essentially straight and parallel. Effective 
separation of the supercritical core and the envelope continues for At > 0, as can be seen 
by the dramatic drop (by more than four orders of magnitude) in M for r ^ _R cr it(At,) (see 
Fig. 6g). As noted in § 1, Safier et al. (1997) modeled the evolution of a cloud's envelope 
during the post-PMF epoch. Equation (104) of Safier et al. gives their predicted accretion 
rate as a function of the initial cloud mass, radius, and the ratio of the initial mean and 
central density. Evaluating their expression for the parameters used in our typical model 
(M d = 98.3 Mq, Ro = 4.29 pc, and (n n ) /n ntC0 = 0.144) yields a spatially uniform envelope 
accretion rate M = 3.7 M Q Myr -1 , which is close to our calculated accretion rate at the 
core boundary (see Fig. 6g). However, for r ^ 0.26i? — 0.11 pc their predicted value of 
M is several orders of magnitude higher than the one we have obtained in the envelope 
of our typical model. The difference can be attributed to the fact that Safier et al. chose 
to study model clouds for which ionization by the external interstellar UV radiation field 
was unimportant, resulting in a lower degree of ionization and a correspondingly higher 
ambipolar diffusion rate in the cloud envelope. 

Finally, the effect that the tidal gravitational field of the central protostar (see eq. 0) 
has on the vertical structure of the cloud can be seen in Figure 6i. For times At > the 
inner region of the core near the central mass is significantly compressed by the tidal force, 
and Z/r becomes <C 1. Further compression will be provided by the squeezing effect of 
the radial magnetic field component (see eq. ||), particularly at late times in the vicinity 
of ri nil e r . All in all, the disk remains geometrically thin at all radii r > ri nnor for the entire 
duration of the simulation. 



4. Discussion 

4.1. Observational Comparisons and Predictions 

We may compare our result for the protostellar accretion rate during the PMF epoch 
with observations of star-forming cores. As shown in § 3.3, the accretion rate rises rapidly 
early on to M cont ~ 9.4 M Myr -1 for At ^ 10 3 yr (see Fig. 2b) and stays at this value 
up to the formation of the hydromagnetic shock. For At ^ 4 x 10 3 yr the shock is able to 
decelerate the infalling matter, and the accretion rate decreases to M C ent — 5.6 M Q Myr -1 
by At ~ 1.5 x 10 5 yr (the time when M ccnt = 1 M ). Therefore M cen t decreases with 
increasing central mass (see Fig. 2c). This is consistent with estimates of ages and accretion 
rates (oc t~^ e ) for young stellar objects, as deduced from evolutionary diagrams inferred 
from observations of Class and Class I objects (e.g., Saraceno et al. 1996). In particular, 
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the lifetimes of Class objects were estimated in this way to be an order of magnitude 
shorter than those of Class I objects, providing evidence for a decrease in the protostellar 
accretion rate as an object evolves from a Class source to a Class I source (e.g., Andre 
1995; Ward-Thompson 1996). Another argument for a time- dependent accretion rate 
was given by Bontemps et al. (1996), who analyzed the observed CO momentum flux of 
several young stellar objects and found a noticeable decline in the CO flux with decreasing 
circumstellar envelope mass. They suggested that this is indicative of a decrease in the 
protostellar accretion rate (which they assumed to be proportional to the mass outflow 
rate) as an object evolves from Class to Class I. 

It is also of interest to note the 13 CO(J =1 — 0) observations of infalling disks for 
the protostellar candidates HL Tauri (Hayashi, Ohashi, & Miyama 1993) and L1551-IRS5 
(Ohashi et al. 1996). HL Tauri has a mass ~ 0.6 M and a surrounding disk with radius 
~ 1400 AU and mass ~ 0.03 M Q . From the observed kinematics Hayashi et al. derive 
an accretion rate ~ 9 M Q Myr -1 at r ~ 700 AU. The embedded protostar L1551-IRS5 
has a mass ~ 0.5 M and is surrounded by a disk with radius ~ 700 AU and mass in the 
range 3.9 x 10 -2 — 8.1 x 10~ 2 M & . Ohashi et al. deduce an accretion rate in the range 
13 - 26 M Myr" 1 at t ~ 600 AU. These values are comparable to our model results for 
At 5 ^ At ^ At 6 (corresponding to 0.2 M ^ M ccnt ^ 1 M ; see Fig. 6a). During this 
period, 6 M Myr" 1 M £ 9 M Myr -1 for r ^ 500 AU (see Fig. Qg). [Note, however, 
that the temperatures of HL Tauri and L1551-IRS5 are in the range 15 — 50 K, which is 
greater than our assumed value of 10 K and should lead to higher values of M and M cent (t); 
e.g., Shu et al. 1987. For a discussion of how quantities scale with temperature in our 
models, see Basu & Mouschovias 1994.] In Figure 7 we show the mass (M — M ccnt ) of the 
gas surrounding the central point mass in our typical model as a function of r/R for the 
same seven times At,- as in Figure 6. (Taken together, Figs. 6a and 7 may be taken to 
represent the evolution of a protostar from a Class to a Class I object.) For times ^ At 5 , 
the surrounding disk mass spans the range 0.01 — 0.1 M for r ^ 500 AU, which agrees with 
the 13 CO disk masses of HL Tauri and L1551-IRS5 cited above. Finally, we note that the 
age of the oldest part of the molecular outflow from L1551-IRS5 is estimated to be ~ 10 5 yr 
(Bachiller, Tafalla, & Cernicharo 1994). This age is consistent with the time At needed for 
the central mass in our typical model to become ^ 0.3 M (see Fig. 2e). 

Another observational consequence of our model is the magnetic field structure in the 
core after PMF. Figures 66 and Qd show that B rZ ~ B Zfi(l inside the core for the radius 
range 2 x 10~ 4 ^ r/R ^ 4 x 10~ 3 . Hence, there is significant curvature of field lines 
(though, as discussed in § 3.3, there is less bending than there would be if the field had 
remained frozen into the neutrals), with bending angles 9b ~ arctan(_B ri ^/_B ZiCq ) in the 
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range 20° - 50° for 180 AU ^ r ^ 3.5 x 10 3 AU. This type of field geometry may be 
described roughly as having an hourglass shape. In fact, sub-mm polarimetry of the cores 
of W3 IRS5 (Greaves, Murray, & Holland 1994), Mon R2 (Greaves, Holland, & Murray 
1995), and OMC-1 (Schleuning 1998) find field geometries suggestive of an hourglass 
shape on sub-parsec scales. In contrast, field lines remain essentially straight and parallel 
in the magnetically supported envelope (r > 0.1 pc), even after PMF. This agrees with 
polarimetric observations of molecular clouds that indicate well-ordered fields on these 
scales (e.g., Hildebrand, Dragovan, & Novak 1984; Hildebrand 1989, 1996; Novak et al. 
1989; Kane et al. 1993; Hildebrand et al. 1995). 

A unique prediction of our model is the large ion-neutral drift speed that occurs 
during the post-PMF epoch. As shown in Figure 6/, effective ambipolar diffusion following 
PMF yields v D ~ \v n \ ^> C in the inner regions of the core. In our typical model we find 
v D ^ 1 km s _1 for At ^ 2 x 10 4 yr on scales r^2x 10 _4 i? — 180 AU. Large drift 
speeds between neutrals and ions (such as HCO + , HCN + , DCO + , to name but a few) on 
these scales are therefore expected in our model. Detection of such drift speeds (through 
high-resolution observations of HCN or HCO + , say) could be used to observationally 
confirm our model results and to distinguish them from those of nonmagnetic collapse 
calculations (e.g., Shu 1977, Hunter 1977, Foster & Chevalier 1993) or of magnetic collapse 
models that do not account for the effect of ambipolar diffusion (e.g., Tomisaka 1996; Li & 
Shu 1997). 



4.2. Features of the Hydromagnetic Shock 

The properties of hydromagnetic shocks in partially ionized gases have been developed 
extensively by many other authors (e.g., Mullan 1971; Draine 1980; Chernoff 1987; Roberge 
& Draine 1990; Draine & McKee 1993; Smith & Mac Low 1997). Because the ion Alfven 
speed t>A,i = -S Zj0q /(47rm i n i ) 1 / 2 = {m n /m- 1 ) l ^ 2 {n ri /n i ) 1 /' 2 vx is much larger than va, \v n \, and 
\vi\ in our model, we expect the outward-propagating shock that develops after PMF to 
have a magnetic precursor. This is indeed what we find in our simulation: the jump in the 
ion speed V\ and the magnetic field strength B zm typically occurs at a distance of 1 to 3 
computational mesh spacings further away from the symmetry axis than the jump in the 
neutral infall speed v n and the column density a n . The displacement between the locations 



10 In agreement with the results of CM94 and CM95, we find that the magnetic tension force is generally not 
negligible in comparison with the magnetic pressure force at any radius (in either the core or the envelope) , 
both before and after PMF. 
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of the head of the disturbance in the neutral fluid and in the plasma and magnetic field 
decreases at later times. Examination of Figures 56, 5c, and 6e reveals that, in the reference 
frame of the shock, the preshock infall speeds are supersonic, and the postshock speeds are 
also supersonic or just slightly subsonic. Similarly, the ion infall speeds (in the frame of 
the shock) are much less than the ion Alfven speed. Therefore the shock we observe in our 
model is probably best classified as being C-type. Making the approximation that in the 
vicinity of the shock the predominant magnetic stress is that due to the magnetic pressure 
gradient, the ion force equation becomes 



An dr 



(16) 



(see eqs. [28c] and [51] in CM93, which contain additional terms, involving in particular the 
magnetic tension force, that could be used to refine the following simple estimate). This 
yields an approximate shock width 



A 



shk 



B 2 t -7 

z,eq,u 'm-^ 
47TCT n v D 

7.8 x 10 13 



B 



z,eq,d 



B 



- 1 



2,eq,u 



Cr n -, 



2(vJC) 



B 



2,eq,d 



B 



1 



2,cq,u / 



(T/10 K) 1/2 (v A JC) 



m n /2.33 a.m.u.) 1/2 {m/0.1 cm" 3 ) (v D /C) 



1 + 0.067 



B 



z,eq,d 



B 



z,eq,u , 



cm 



(17a) 

(m K2 /2 a.m.u.) 
(mj/30 a.m.u.) 

(17b) 



where -B 2ieq ,u and -B 2ieq ,d are the values of B z e<i upstream and downstream of the shock, 



and v a, u is the upstream Alfven speed. In deriving the last equality of equation (|17a|) 
we have used the relation cr n = 2p n Z; equation (|l|) and the relation C = (k B T ' /m n ) 1 ^ 2 
have been used in deriving equation ( 17b ). At the time Ate the shock front is located 



at r ~ 3.9 x 10 3 Rq = 5.2 x 10 16 cm, and, in the vicinity of the front, n\ ~ 



10 2 cm 3 , 



t> A u ~ C, B z ^ d / B z ^ n ~ 4.6 (see Fig. 66), and v n ~ 0.4C (see Fig. 6/). For these values 
our rough estimate for the shock width given by equation ( |17b| ) yields A s hk — 4.3 x 10 16 cm. 
Examination of Figure 6e at the time A^6 reveals that the shock has an actual width 
A sh k — 2.1 x 10~ 3 i? = 2.8 x 10 16 cm. (Thus A shk /r sh k ~ 0.5 at that time, so the thin-shock 
approximation that underlies the estimate [17] is marginally satisfied.) 



n A C-type shock is characterized by neutral velocities that (in the shock frame) remain supersonic 
throughout. Hence the shock in our simulation cannot be strictly of this type when the downstream neutral 
speed is subsonic. In that case a viscous (J-type) subshock may form (Draine & McKee 1993), although, as 
noted by Li & McKee (1996), turbulent diffusivity behind a real shock could plausibly keep the postshock 
flow supersonic and thereby obviate the need for such a subshock. 
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As noted in § 1, Li & McKee (1996) proposed that a hydromagnetic shock would form 
in a collapsing core as a result of the decoupling of flux from the ion and neutral fluids 
because of Ohmic dissipation (a process that becomes important in regions of density 
n n ^> 10 11 cm~ 3 ; see § 2.1). They argued that the accumulating flux diffuses outward to 
regions of lower density, where improved coupling with the matter causes it to present an 
obstacle to the infalling neutral gas — thereby giving rise to a hydromagnetic shock. While 
our numerical results have confirmed Li & McKee's basic shock-formation scenario, they 
have revealed that ambipolar diffusion, which mediates the shock, can also, following PMF, 
halt the inward advection of magnetic flux on scales (r ^ 5 AU) where Ohmic dissipation 
is not important. In other words, our results have shown that, during the post-PMF 
epoch, the field-matter decoupling that drives the hydromagnetic shock is due to ambipolar 
diffusion alone and does not depend on the effect of Ohmic dissipation at r < r inner . 

Because of the similarity in the basic shock-formation mechanism, it is of interest 
to compare our detailed simulation results with the predictions of the (simplified and 
analytic) shock model of Li & McKee (1996). From their requirement that the magnetic 
pressure of the shock balance the ram pressure of the neutrals, which were assumed to be 
freely-falling into the shock, Li & McKee derived relations for the shocked magnetic field 
strength and the shock location (see their eqs. [7] and [8]) in terms of the accretion rate M, 
the flux-to-mass ratio in units of the critical value for collapse (dubbed e in their paper), a 
parameter related to the logarithmic gradient of the magnetic field (dubbed x), the ratio 
Z/r of the local gas scale height and the radius (dubbed h), and the protostellar mass 
(denoted by m*). At the time At 6 we have at the location of our shock M 9 M Q Myr -1 , 
e 0.9, x ~ 2, h ~ 0.3, and m* = M ccnt (/S.t§) fa 1 M . Inserting these values into their 
equations (7) and (8) yields a shocked magnetic field strength ~ 630 fiG and a shock radius 
~ 2.1 x 10 3 AU; by comparison, in our model the shocked magnetic field strength at that 
time is -B 2j eq,d ~ 330 fiG and the shock radius is r s hk ~ 3.5 x 10 3 AU (see Figs. 6b and 
6e). The main reason why the analytic expression overestimates the numerically calculated 
magnetic field strength is that, contrary to the assumption of Li & McKee, the preshock 
neutrals are not in free fall but, rather, are strongly decelerated by magnetic forces (in our 
simulation we find that the preshock acceleration of the neutrals is reduced to ~ 0.25g r ). 
The corresponding reduction in the preshock ram pressure leads to a lower value of the 
postshock field amplitude, with a further reduction in the calculated field strength brought 
about by the contribution of magnetic tension (ignored in the analytic estimate) to the total 
magnetic force. Since the analytic estimate of r s hk is based on relating the postshock field 
strength to the total magnetic flux inside the shock, the overestimate of the field strength 
naturally results in an underestimate of the shock radius. Despite these discrepancies, Li & 
McKee's analytic representation of the shock parameters provides a decent approximation 
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to the results of our numerical calculation. 

Comparison of the column density a n (see Fig. 6c) and neutral infall speed v n (see Fig. 
6e) behind the shock shows qualitative agreement with Figures 2a and 2d of Li Sz McKee 
(1996), including the free-fall behavior near the central protostar. In the pre-shock region, 
however, our infall speeds are smaller than theirs because of their assumption of free fall 
upstream of the shock: Li & McKee typically overestimate v a>u by a factor ~ 2. As a result, 
the Mach number of the shock relative to the upstream gas (and thus the shock strength) is 
greater in their model than in ours (see Fig. 5c) by a similar factor. We have not compared 
our results for the magnetic field structure behind the shock with those of Li & McKee on 
account of the fact that their system of MHD equations was not closed (it did not include 
the magnetic induction equation), so that they were unable to calculate the magnetic field 
with any accuracy (see their Fig. 2b). 



4.3. Stability of the Core Against Magnetic Interchange 

Our simulation has revealed that rapid ambipolar diffusion occurs behind the outward- 
propagating HMD. The effect that this has on the mass in the flux tubes downstream of the 
HMD can be seen in Figure 8a, which shows the local mass-to-flux ratio dM/d^B — &n/B Zjeq 
(normalized to the critical value for collapse) as a function of t/Rq for the same seven times 
Atj as in Figure 6. For At > At\ a local minimum in a n /B z ^ eq appears after the passage of 
the HMD. Hence, there is a region behind the HMD for which d(a n /B Z)eq )/dr > 0. This is 
a necessary condition for the onset of an interchange instability (e.g., Spruit & Taam 1990; 
Lubow & Spruit 1995; Spruit, Stehle, & Papaloizou 1995). Li & McKee (1996) speculated 
that such a situation could arise in the wake of a hydromagnetic shock in a collapsing core 
and suggested that it would act as source of turbulence in the postshock region of the 
inflow. 

Blaes & Balbus (1994) found that the magnetic shearing instability in differentially 
rotating disks could be stabilized if the disk is weakly ionized. This will also be the case 
for the interchange instability in a weakly ionized disk: instability is possible only if the 
growth rate 7 n and the neutral-ion collision time r ni satisfy the condition 7 n r ni < 1. This 
condition reflects the fact that there has to be sufficient collisional coupling between the ion 



12 Li & McKee (1996) also noted that the shock may be unstable to the Wardle instability, which involves 
ions collecting in magnetic field "valleys" (Wardle 1990). However, the shock will be immune to this 
instability if the ion density is determined by the local chemical reaction balance (as assumed in our 
calculation) rather than by the divergence of the ion mass flux. 
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and neutral fluids for a magnetic interchange instability to grow in the neutrals; otherwise 
the instability is damped. Calculation of 7 n in our model is hampered by the fact that 
previous studies of interchange instability have been carried out only for disks that are in 
hydrostatic equilibrium, with exact balance between magnetic and gravitational forces. We 
can apply the results of these studies to our model only if the region behind the shock where 
d(a n /B z eq )/dr > is effectively in quasi equilibrium, with approximate balance between 
gravitational and magnetic forces, and with infall speeds |t> n | <C {r\g r \) 1 ^ 2 (~ the free-fall 
speed). In our model, the magnitude of the acceleration of the neutrals in this region of the 
core does not become ^ 0.1|<jv| until times ~ At 6 ; hence, approximate equilibrium between 
magnetic and gravitational forces is valid only at these later times. Spruit & Taam (1990) 
found that the growth rate for the most unstable linear interchange modes is 



The product 7 n T n j for the region of the core susceptible to interchange instability is shown 
in Figure 8b as a function of t/Rq at the time At 6 . We also plot (Fig. 8c) the product 
Tii^kin, where Tki n = r/\v n \ is the kinematical timescale, as a function of r/i? for the same 
region of the core and time. If this product is < 1, the unstable modes will be "swept up" 
by the infalling gas before they have time to grow. It is evident from these figures that 
7 n T ni < 1 and 7 n r kin > 1 for the region of the core susceptible to magnetic interchange. This 
means that there is sufficient collisional coupling between the ions and neutrals, and that 
the instability will grow before being swept along with the neutrals. Hence, this region of 
the core may be interchange unstable. An instability of this type would enhance the tansfer 
of gas with a high mass-to- flux ratio to the center (e.g., Spruit & Taam 1990), and, as 
noted by Li & McKee (1996), might also lead to the development of turbulence that could 
increase the field diffusivity in the postshock gas. However, the onset and development of 
this instability can only be studied by means of a fully 3-D simulation. 



The magnetic flux problem in star formation has to do with the fact that the magnetic 
flux of a 1M blob of matter in the diffuse interstellar medium is typically several orders of 
magnitude greater than the flux of a 1M protostar. Such a blob of matter would therefore 
have to get rid of most of its flux before becoming a star. Ambipolar diffusion has long 
been suggested as a means by which the magnetic flux problem could be resolved (e.g., 
Mestel & Spitzer 1956; Mouschovias 1978; Paleologou & Mouschovias 1983; Nakano 1984; 
Mouschovias, Paleologou, & Fiedler 1985). In general, these earlier studies focused primarily 
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on the role of ambipolar diffusion and the magnetic flux problem for the pre-PMF phase of 
protostellar evolution. P] While ambipolar diffusion does indeed reduce the flux-to-mass 
ratio during that phase, the flux contained within a 1M & region of a molecular cloud core 
is still much larger at the time of PMF than typical protostellar fluxes. Specifically, CM94 
and CM95 found that the central 1M flux tube within their cores had a total magnetic 
flux ~ 10 30 Mx (consistent with our results in § 3.2) during the pre-PMF dynamical collapse 
phase of the typical model. This value represents a reduction by a factor ~ 5.6 of the flux 
associated with that mass before the onset of ambipolar diffusion. Nevertheless, it greatly 
exceeds the plausible upper limit (~ 6 x 10 26 Mx) on the flux of a solar-mass protostar 
(estimated assuming an average surface field of 10 kG and a stellar radius of 10 11 cm; see 
Li & McKee 1996). 

We have shown in this paper that the rate of ambipolar diffusion is strongly increased 
during the post-PMF epoch of star formation. It is therefore of interest to examine the 
implications of our simulation results to the magnetic flux problem. As discussed in § 2.2, 
in our calculations we only consider the core regions at radii r ^ r inner (~ 7.3 AU for our 
typical model), where ambipolar diffusion is the dominant mechanism of flux loss. Initially, 
the magnetic flux contained within r inner is ^B.cent {At = 0) = 6.7 x 10 26 Mx. As shown in 
Figure 3a, the central flux increases before the onset of rapid ambipolar diffusion. This 
continues to the time At ~ 10 3 yr. For At > 10 3 yr, ambipolar diffusion prevents further 
advection of flux from r > r innor into the central sink, and $B iCen t changes very little after 
this time. By the time At « 10 5 yr, when M cent « 1M , $£ jCcnt ~ 5 x 10 27 Mx. This 
represents a decrease of over two orders of magnitude relative to the flux associated with 
this mass at the time of PMF. While this value is still about an order of magnitude higher 
than our adopted upper limit on the protostelar flux, the discrepancy is now much lower 
than previous estimates of ambipolar diffusion have indicated. The important conclusion 
from our work is thus that ambipolar diffusion in contracting molecular cloud cores can in 
principle contribute significantly to the resolution of the magnetic flux problem by reducing 
the magnetic flux brought into a solar-mass protostar by a factor ^ 10 3 . The new, and 
somewhat surprising, result is that most of this reduction occurs after PMF. 



13 On the basis of a consideration of the timescales for ambipolar diffusion and Ohmic dissipation at high 
densities, Nakano & Umebayashi (1986b) suggested that significant flux loss could only take place (primarily 
by Ohmic dissipation, according to their estimates) during the dynamical phase of core collapse. Lizano & Shu 
(1989) similarly concluded that the resolution of the protostellar magnetic flux problem must occur during the 
dynamical stage of core evolution: using the quasi-static approximation (valid for n n c <^ a few x 10 cm ; 
Fiedler & Mouschovias 1993; CM94; Basu & Mouschovias 1994) to calculate the contraction of a slightly 
subcritical molecular cloud, they found that only a small amount of flux is lost by ambipolar diffusion from 
the central flux tubes before runaway collapse is initiated. 
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The remainder of the protostellar magnetic flux could possibly be extracted from the 
infalling mass through Ohmic dissipation within r inner , although refreezing of the magnetic 
field into the matter, brought about by collisional reionization at densities n n ;> 10 14 cm~ 3 
(e.g., Pneuman & Mitchell 1965; Nakano & Umebayashi 1986b; Li & McKee 1996), as 
well as anomalous diffusivity (e.g., Norman & Heyvaerts 1985) operating in the reionized 
gas, could complicate the issue. Another complicating factor is the strong likelihood that 
much of the mass and flux carried into the protostar pass through a rotationally supported, 
circumstellar accretion disk of size 3> r inncr (e.g., Lubow, Papaloizou, & Pringle 1994; 
Reyes- Ruiz & Stepinski 1996; Li 1996). It is also conceivable that magnetic flux is brought 
to the vicinity of the protostar but excluded from its interior by turbulent diffusivity 
associated with convection. Since the region within r inner was excluded from our calculation, 
we do not pursue this topic any further in this paper. 

5. Summary 

We have simulated the formation and growth of a central (i.e., protostellar) point 
mass in a gravitationally collapsing core of a nonrotating magnetic cloud modeled as an 
isothermal and axisymmetric thin disk. Following up on the results of previous simulations, 
we concentrated in this study on the core evolution after point-mass formation (PMF), 
paying particular attention to the role of ambipolar diffusion during this phase. In view of 
our assumptions of gas isothermality and magnetic flux freezing into the ions, our results 
are applicable only on scales r ^ 5 AU. 

Just prior to the formation of a point mass, the model core is dynamically collapsing 
(though not freely falling), with infall speeds that are comparable to, or exceed, the 
isothermal speed of sound and the local fast-magnetosonic speed, and accelerations ranging 
from 0.25 to 0.5 of the local gravitational acceleration. 

We have calculated the evolution of the central protostar up to the time that it has 
grown in mass to 1M Q . Ambipolar diffusion causes the evolution of our model core to differ 
significantly from that found in previous calculations of dynamically collapsing cores, which 
considered either nonmagnetic or magnetic but perfectly conducting clouds. In particular, 
we find that ambipolar diffusion in the weakly ionized gas surrounding the central protostar 
is "revitalized" by the increase in the strength of the gravitational field brought about 
by the formation and growth of the central point mass. (An alternative, but equivalent, 
explanation of this revitalization is that the magnetic tension force acting on the ions 
increases to the point where the ion-neutral drift speed becomes comparable to the neutral 
inflow speed.) The ambipolar diffusion becomes rapid enough to stop the inward advection 
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of magnetic flux, which therefore begins to pile up in the inner regions of the collapsing 
core. The region of piled-up flux develops into a hydromagnetic disturbance that propagates 
outward and increases the local magnetic field amplitude away from the axis of symmetry. 
As the field increases and the disturbance reaches regions where the ion-neutral coupling 
becomes stronger, the hydromagnetic front evolves into a (C-type) shock: for the physical 
parameters assumed in our typical model, this occurs at radii r ^ 90 AU. The shock speed 
is supersonic (and super fast-magnetosonic) with respect to the infalling upstream gas, but 
for r ^ 140 AU it becomes subsonic in the protostellar frame. The shock decelerates the 
neutrals and interrupts their infall, thereby decreasing the accretion rate onto the central 
protostar: we confirmed that the accretion rate obtained by accounting for the effect 
of ambipolar diffusion was less than that calculated in a model that evolved under the 
assumption of flux freezing into the neutrals. Far behind the shock (at radii ~ 10 — 40 AU 
in our typical model) the neutrals are reaccelerated to free-fall, and the column density and 
neutral infall speed scale as r" 1 / 2 . The column density scales with the time At since PMF 
as cr n oc (At) s At, with s At lying in the range —0.76 £ s At ^ —0.55 (whose upper bound 
is close to the value s At = —0.5 obtained under the assumption of self-similarity). Our 
numerical results are in overall agreement with the simplified analytic model of Li & McKee 
(1996), who first predicted the existence of the outward-propagating C-shock. We have 
found that, after PMF, this shock is driven primarily by ambipolar diffusion, which leads 
to field-matter decoupling on scales larger than those where Ohmic diffusivity effects are 
important. CCK97 derived a semianalytic similarity solution that incorporates ambipolar 
diffusion and captures the main features of the post-PMF core evolution (including the 
shock formation) found in our simulation. 

For the typical model presented in this paper, the protostellar accretion rate increases 
from ~ 5 Mq Myr" 1 just prior to protostar formation to a maximum ~ 9.4 M Myr -1 
at At ~ 10 3 yr (although, because of the approximations involved in our calculation of 
the radial magnetic field component at the inner boundary, the actual maximum accretion 
rate might be higher). The accretion rate subsequently decreases, largely on account of the 
interruption of the infall by the hydromagnetic shock, and it is equal to 5.6 M Myr -1 by 
the time (1.5 x 10 5 yr after PMF) the central mass has grown to 1 M . For comparison, the 
"canonical" accretion rate (~ C 3 /G; Shu 1977) for this model is 1.5 M Myr" 1 , although 
it has been recognized that this value would be larger in clouds where magnetic stresses 
supplement thermal pressure support. [For example, Li & Shu (1997) estimated an increase 
by a factor ~ (1 + H ), where H (~ 1) is the magnetically supported fractional overdensity 
in the stationary, pre-PMF cloud. Other magnetic models (e.g., Tomisaka 1996; Basu 1997) 
have also indicated a larger-than-canonical average accretion rate.] The decline in the 
accretion rate with time is consistent with the trend inferred from observations of Class 
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and Class I sources. Furthermore, the accretion rate and the circumstellar mass at radii 
^ 500 AU in our typical model are comparable to those derived by millimeter observations 
of protostars with "infalling" disks, such as HL Tauri and L1551-IRS5. The magnetic field 
structure in our model is consistent with polarimetric observations of molecular clouds 
that reveal essentially uniform fields in cloud envelopes (r ^ 0.1 pc), and with sub-mm 
polarimetric surveys that indicate hourglass field structures within cloud cores. 

Ambipolar diffusion is so effective after PMF that it leads to a reduction by more than 
two orders of magnitude in the flux threading the central 1 M Q , compared with less than 
one order of magnitude reduction in the flux threading this mass between the start of the 
cloud contraction and the time of PMF. Altogether, ambipolar diffusion reduces the flux by 
more than three orders of magnitude and brings it to within an order of magnitude of the 
estimated upper limit on the flux of a ~ 1 M & protostar. Ambipolar diffusion in collapsing 
cloud cores could thus go a long way toward resolving the protostellar magnetic flux 
problem, although to fully address this issue one would need to consider Ohmic diffusivity 
on scales smaller than those included in the present calculation as well as the effects of 
rotation. 

Large 1 km s _1 ) ion-neutral drift speeds occur in our representative simulation 
on scales r <^ 200 AU for At ^ 2 x 10 4 yr, and are a unique prediction of our model. 
Observational detection of this effect could be used to distinguish our results from other 
(nonmagnetic and magnetic) collapse models that do not account for ambipolar diffusion. 
We have also confirmed that the postshock region is susceptible to the magnetic interchange 
instability, as first suggested by Li & McKee (1996). This instability could only develop 
outside the region of large ion-neutral drift velocities, where it is not subject to quenching 
by ambipolar diffusion effects. Magnetic interchange would increase the rate of mass transfer 
to the center and could also lead to turbulence that might enhance the magnetic diffusivity. 
A numerical investigation of this instability does, however, require a 3-D simulation. 

The mass and magnetic flux distributions during the post-PMF epoch are directly 
relevant to the question of protostellar disk formation and the generation of disk outflows 
by magnetic ejection mechanisms (e.g., Konigl & Ruden 1993). However, in order to model 
the formation of centrifugally supported circumstellar disks, it is necessary to incorporate 
rotation and magnetic braking into the core-evolution calculations. This will be considered 
in a future publication. 
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APPENDIX A 

DESCRIPTION OF MODEL PARAMETERS 

We consider isothermal, magnetic interstellar molecular clouds consisting of neutral 
molecules (H 2 with a 20% He abundance), singly charged molecular and atomic ions (such 
as HCO + , Mg + , Na + , Fe + , C + , S + , and Si + ), and electrons. For the purposes of an accurate 
calculation of the equilibrium abundances of charged particles, we also include (negative) 
singly charged grains and neutral grains. However, the collisional effects of grains on the 
neutrals (which in certain cases can be significant; see CM94, CM95) are ignored in this 
paper: this is done by setting the parameter {aw) gn (see below) equal to 0. The effects 
of magnetic braking and rotation, which can also affect the evolution of a core (Basu & 
Mouschovias 1994, 1995a, b) are similarly neglected in this paper. They will be accounted 
for in a later publication. 

The assumptions that model clouds are thin and that evolution along field lines is 
quasistatic allow the time-dependent nonlinear set of nonideal MHD equations that govern 
the evolution of the model cloud to be integrated over z, thus reducing the dimensionality 
of the formal problem (CM93). The simplified set of equations is listed as equations 
(Al)-(A17) in CM95. The equations are cast in dimensionless form by adopting the 
quantity 2iiGa c ^ e f as the unit of acceleration, B Te { as the magnetic field strength, and the 
isothermal speed of sound C as the unit of velocity. The quantity <7 C)re f is the central column 
density of a reference state used to specify how magnetic field lines are initially loaded 
with mass, and B Tc{ , as mentioned in § 2.2, is the background magnetic field strength at 
infinity. The implied units of time, length, and mass density are, respectively, C/2nGa c ^ e f, 
C 2 /2-kGo- C:IC {, and 2^Gal rei /C 2 . The resulting system of dimensionless equations contains 
the following set of nondimensional parameters: 

{dM/d^BU KrcfMcf) _ n 1 Qfi / <r C ,rcf \ / 30 (J.G \ 

^ " (dM/d* B ) A:C , lt ~ (,,_>:; v'77) " 1.3.63 x 10-3 g cm^j I - J 



5 ^cxt _ / P cx t \ /3.63 x 10-3 g cm" . 

= (tt/2) Go* ~ ° A Vl.38 x 10"" dyn cm"' J [ ^~ »" {A2) 
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where the quantity a in equation ( |A4| ) is the grain radius. The parameter /id.co is the initial 



central mass-to-flux ratio in units of the critical value for collapse. P ext is the constant 
external pressure normalized to the self-gravitational stress of the matter contained within 
the central flux tube of the reference state. The parameter (aw)n{ 2 is the dimensionless 
neutral-ion collisional rate, and {aw) gn is the dimensionless neutral-grain collisional rate. 
The initial equilibrium state introduces one more parameter, the quantity Z re f, a length 
scale in the column density of the reference state (see CM93, eqs. [60b] and [73b]). The 
dimensionless cloud radius is taken to be 5 Z re f in this paper. A final dimensionless parameter 
is the initial dust-to-gas mass ratio, Xg,o ; which is equal to 0.01 in the typical model. 

Abundances of charged particles are determined by solution of the chemical rate 
equations with balance between creation and destruction of the various charged species 
specified above, accounting for such processes as cosmic-ray ionization, dissociative 
recombination of molecular ions and electrons, radiative recombination of atomic ions and 
electrons, attachment of charged particles on grains, and charge transfer between molecular 
and atomic ions; we also account for ionization due to an external (interstellar) ultraviolet 
radiation field. The relevant chemical equations are given by equations (7)-(17) of CM95. 
They contain six relevant nondimensional parameters of the form C oU vcr = C oU v/C C r' 
where ( ao is the UV ionization rate at the cloud boundary for neutral species a , and C CR is 
the cosmic-ray ionization rate. 

APPENDIX B 



DENSITY SCALING IN THE INNER FLUX TUBES OF MODEL CLOUDS 

For our thin-disk model, the density p n (= m n n n ) is calculated from equation @. We 
examine the power-law behavior of the density in the inner flux tubes of the core for the 
case in which the tidal gravitational field of the central protostar is negligible, and then for 
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the opposite case where it is the dominant term in equation (0). In the first case, using 
the fact that P cxt <C (7i/2)Ga 2 in the inner flux tubes of the core and that the magnetic 
squeezing term is also negligible there for most times of interest (see discussion in §2.2 
following eq. 0), the above equation yields p n = 7iGa 2 /2C 2 . During the later stages of the 
pre-PMF phase of core contraction, cr n oc r _1 (CM94, CM95; Basu & Mouschovias 1994, 
1995a,b; see also § 3.2). Therefore p n oc r~ 2 as the core approaches PMF. 

We now turn to the case where the tidal stress from the central protostar is the 
dominant term on the r.h.s. of equation (|9|). Near the protostar, tidal squeezing of the 
disk (see Fig. 6i) results in Z/r = a n /2p n r < 1. Q Expansion of the second term in 
braces on the r.h.s. of equation (0) gives the solution p n = (GM cent /8) 1/2 a n /CV 3 / 2 . In the 
limit of free-fall collapse near the protostar a n oc r -1 / 2 (see § 3.3) and p n is again oc r~ 2 . 
Hence, the power-law behavior of the density in our model is constrained to be oc r -2 
before and after PMF. This result is different from the behavior p n oc r~ 3 / 2 , which would 
occur for the case of spherical infall in the gravitational field of a central point mass. This 
difference is due to our use of the vertical one-zone approximation for the mass density 
[p n (z,r,t) = p n (0,r, t) = p n (r, t)} in equation (|9|) for balance of forces along field lines. 
Figure 9 shows the number density n n , normalized to n n)C o = 2.6 x 10 3 cm -3 , as a function 
of r/Ro at the same seven times At,- as in Figure 6. The r~ 2 scaling of the density is 
apparent in this figure at small radii. 

APPENDIX C 

BREAKDOWN OF FLUX FREEZING AFTER POINT-MASS FORMATION 

In this Appendix we present a heuristic argument that demonstrates the breakdown 
of flux freezing, and the corresponding revitalization of ambipolar diffusion, following 
point-mass formation (PMF). We show this by means of a reductio ad absurdum argument: 
we assume that the collapse proceeds under flux freezing conditions, and we use the 
frozen-flux solution to evaluate the ratio of the ion-neutral drift speed (v D ) and the free-fall 
speed (v B = \2GM{r) /r} 1 / 2 ). Flux freezing corresponds to v D <C v B ; however, by deriving 
the magnetic force on the ions from the numerically computed evolution and using the 
ion equation of motion to calculate the resulting ion-neutral drift speed, we find that the 
ratio v D /v g becomes 3> 1 after PMF. This implies that the assumption of flux freezing 
is not self- consistent and that ambipolar diffusion necessarily sets in. According to this 



14 This squeezing will be augmented, particularly at late times, by magnetic field effects represented by the 
second term on the r.h.s. of equation (H). 
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argument, the increase of v D /v s in the inner core following PMF can be attributed to the 
growth of the magnetic tension term in the ion force equation. The quantity v D /v B is 
essentially the inverse of the ratio of the ambipolar diffusion timescale (r AD = r/v D ) and 
gravitational contraction timescale [r gr = (r/\g r \) l l 2 « r s ] considered in § 3.3 (see eq. |T5fl). 
In the argument presented in § 3.3 (see also § 1), T AD /r gr is shown to decrease below 1 
during that epoch on account of the effect of the central point mass on the magnitude 
of the gravitational acceleration in its vicinity. As we show below, these two alternative 
descriptions of the ambipolar-diffusion revitalization process are equivalent. 

In a disk-like cloud, the dominant terms in the expression for the radial magnetic force 
per unit mass are 

-^mag.r _ 1 | rj r> y^^eq | 

a n 27ra n V or J 




9Bz,eq 



2na n \ B ZtCq dr 



(C5) 



(see eq. [28c] of Ciolek & Mouschovias 1993). The first term on the r.h.s. of equation 
is the magnetic tension force, while the second term is the magnetic pressure force. It turns 
out (see below) that, in a frozen-flux core, the tension term comes to dominate after PMF. 
In that case, using the ion force equation 



— v D = F magjr (C6) 



• m 



and the continuity equation (M = — 27rra n v n ), and assuming \v n \ ~ v s in the inner flux 
tubes of a collapsing core, one obtains 



Tni M 



, (C7) 



_2fi 2 B M 

where /is is the mass-to-flux ratio in units of the critical value for collapse, as discussed 



in § 1. The term in brackets on the r.h.s. of equation (|C7 ) is the same as equation (4) 
of Li & McKee (1996), who used e = /i^ 1 in their expression. Consider now the behavior 
of B r ^z and B Zfiq in a frozen-flux core following PMF. As was noted in § 2.2, B ZfiCL oc o n 
and B r} z oc g r in a thin, perfectly conducting disk with \ib = const. Thus, outside the 
central sink cell, B r Z ~ $b(?"> t)/2nr 2 w & BtCcnt /2nr 2 after PMF (see eq. and associated 
discussion). Here $#(r, t) is the total flux enclosed within radius r, which is ~ $_B iCen t 
after PMF. Hence, B r z scales as r~ 2 . Near the central point mass the inflowing matter is 
in approximate free fall and hence B ZtBq oc a n oc r" 1 / 2 (see § 3.3). Therefore, after PMF, 
B r Z /B zm oc $s,cent/ r3 ^ 2 will become 3> 1 as r — > 0. This behavior has indeed been seen 
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in our frozen-flux model simulations, and is also apparent in the self-similar frozen-field 
collapse calculations of Li & Shu (1997) and CCK. 

Because of the increase of B ry z/B z>eq near the origin after PMF, the ratio v D /v B given 
by equation (C7) will also become 3> 1 there. This is illustrated graphically in Figures 10a, 
106, and 10c, which show, respectively, the ratios v D /\v n \, B^ z /B Zfiq , and -F m a g ,r/-Fmag,pres in 
the inner flux tubes of the frozen-flux model core presented in § 3.3. The quantity -F mag ,r is 
the actual magnetic force per unit mass calculated in our model (given by eq. ||C5|| above), 
which includes magnetic tension, whereas -F magiPres = B 2 /2na n is the force term used by 
Li & McKee (1996), which approximately represents the effect of magnetic pressure alone. 
It is seen from these figures that v D /\v n \ does indeed become 3> 1 as the core evolves, and 
that this is due to the increase of B r ^ z j B Zfi ^ which scales as r -3 / 2 , as expected. (As pointed 
out in § 3.3, our frozen-flux model is calculated by setting v\ = v n , so that v D = in our 
numerical code; to obtain v D in Fig. 10a, we used equations [|C6|| and ||C5|| .) These results 
also confirm that -F m ag,r/-^ma g ,pres = B r ^z/B Z)C<l . The fact that f D /|f n | becomes ^> 1 means 
that the assumption of continued flux freezing is not self- consistent and that ambipolar 
diffusion must eventually take place. []] The region of flux decoupling will move outward 
with time, which is consistent with the fact that ambipolar diffusion in our simulations 
first becomes noticeable on our smallest resolvable scale, and that the decoupling front 
subsequently moves outward. 

This u B T z/ B ZtBq " explanation of the onset of ambipolar diffusion is the same 
as the "r gr /r n i" argument given in §§ 1 and 3.3, where we emphasized the fact that 
T Ao/ T gr ~ TgrAiu- This can be seen in the following way. Substituting B ri z ~ $s(r, t)/27rr 2 
into the ion force equation (|C6|), we get 

V Ii = -^-B xeQ -^ . (C8) 

2no- n * q 2vrr 2 V ; 

Dividing this equation by v s , and using the relation M = v^r/2G, yields 



47r 2 Ga n M 



2r y/2n 2 B r gr 



r '" (C9) 



where, in the last equality, we used r gr = y/2r/v s and the definition of From this 
equation it follows that T AD /r gr w v a / v v = A*B r gr/ r ni- The expression given by equation 



15 This conclusion did not follow from the estimate of v D /v a in Li & McKee (1996) because they took the 
radial magnetic force per unit mass to be -F m ag,pres rather than F magjr . and thus omitted the ratio -B r ,z/-Bz,eq 
from their equation (4). The argument presented in this Appendix also explains why the revitalization of 
ambipolar diffusion following PMF cannot be studied using a spherically symmetric model (e.g., Li 1998), 
where only the magnetic pressure force, but not the tension force, is taken into account. 
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( P9|) is the same as that derived by Mouschovias (1991, eq. [10a]), who found that 
r AD oc (T| r /r ni )( ( l ) BiCI . it /$ B ) 2 , where ^b,cth is the critical magnetic flux. Finally, by inserting 
equations ( PB| ) and ( |C£| ) into equation ( JUty and using the relations \g r \ ~ GM/r 2 and 
B T z ~ $_B/27rr 2 in the central flux tubes of a core following PMF (see eqs. || and ||), one 
obtains (neglecting thermal-pressure forces) 

(This relation can also be derived from eq. [24] of Basu 1997.) Therefore, equations ( |C9"| ) 
and QUlPp yield equation (|T3p, and the u B r ^ z / B z eq " and "r gr /r ni " pictures are seen to be 
completely equivalent, as claimed. 
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Table 1 

DlMENSIONLESS PARAMETERS OF A TYPICAL MODEL CLOUD 1 " 



A^d.cO 


^ref 


-Text 


(<7w)iH 2 




0.256 


5.57T 


0.1 


2.53 x 10" 48 






t For a discussion of the meaning of these parameters, see Appendix A. 

Table 2 

DlMENSIONLESS UV IONIZATION PARAMETERS tt 



Cc ,UV,CR 


CvIg ,UV,CR 


C-Jao.UV.CR 


Cfc ,uv,cr 


Cs ,UV,CR 


Csi ,UV,CR 


(10 6 ) 


(10 6 ) 


(10 6 ) 


(10 6 ) 


(10 6 ) 


(10 6 ) 


2.60 


0.90 


0.11 


2.40 


14.4 


24.0 



All ionization rates are normalized to a cosmic-ray ionization rate of 5 x 10~ 17 s _1 . All 
dimensional UV ionization rates are taken from Table 9 of Black & Dalgarno (1977, and 
references therein). 
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Figure Captions 

Fig. 1.— Spatial profiles of physical quantities in the typical model, prior to point-mass 
formation, as functions of radius r (normalized to the initial cloud radius Rq 
= 4.29 pc) at eleven different times tj (j = 0, 1, 2, 10 ) chosen such that the central 
density has increased by a factor of lCP with respect to its initial value. These times 
are 0, 7.56 Myr, 9.29 Myr, 9.56 Myr, 9.607 Myr, 9.6167 Myr, 9.6169 Myr, 9.6190 Myr, 
9.61977 Myr, 9.61981 Myr, and 9.61983 Myr, respectively. An asterisk on a curve, 
present only when a supercritical core has formed, marks the instantaneous position 
of the critical flux tube. An open circle on every curve locates the instantaneous 
position of the critical thermal lengthscale. As discussed in § 2.2, the assumption of 
isothermality and the form of the induction equation used in our calculation are not 
valid for densities n n > 10 11 cm -3 m 4 x 10 7 n n)C o. (a) Neutral density, normalized 
to its initial central value n ri)C o (= 2.6 x 10 3 cm -3 ), (b) Neutral column density, 
normalized to its initial central value <7 n;C o (= 5.59 x 10 -3 g cm" 2 ), (c) Mass-to-flux 
ratio, in units of the critical value for gravitational collapse. Between times t\ and 
t 2 a critical core has formed (because of ambipolar diffusion) inside a magnetically 
subcritical envelope, (d) Vertical (z) component of the magnetic field in the equatorial 
plane of the disk, normalized to its initial central value -B 2)eq ,co (= 35.3 /xG). (e) Radial 
(r) component of the magnetic field at the upper surface of the cloud, normalized to 
B Zt eq,co- if) Infall speed of the neutrals, normalized to the isothermal speed of sound 
C (= 0.19 km s -1 ). By the end of the run, |t> n | m C in the inner core, (g) Ion-neutral 
drift speed (= v\ — v Q ), normalized as in (_/). Ion depletion at higher densities is 
responsible for the appearance of the maximum in the core for t > t$. (h) Mass infall 
rate, in units of M Myr -1 . (i) Ratio of cloud vertical half-thickness Z(r) and radius 
r. 

Fig. 2.— Accretion in the central cell during the PMF epoch, (a) Central (protostellar) 
mass, in M , as a function of time At (= t — tio). Formation of a point mass occurs 
at At <; 200 yr. (b) Accretion rate in M & Myr -1 . Dashed curve refers to a model that 
did not include the effect of ambipolar diffusion for At > 0. For At ^ 4 x 10 3 yr, 
ambipolar diffusion, through the formation of a hydromagnetic shock in the neutrals, 
decreases the protostellar accretion rate of the typical model, (c) Accretion rate 
in M Q Myr -1 as a function of central mass (in M & ). Dashed curve refers to the 
frozen-flux model, as in (b). 

Fig. 3.— Physical quantities in two different computational mesh cells, as functions of At 
(pd the time elapsed since PMF). (a) Cell 1: ratio of ambipolar-diffusion timescale 
t ad and gravitational contraction timescale r gr , ratio of gravitational timescale and 
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neutral-ion collision timescale r ni , and (normalized) magnetic flux $ B contained 
within the cell. The gravitational contraction timescale decreases because of an 
increase in the magnitude of the gravitational field associated with the central 
accreting protostar, which increases the rate of ambipolar diffusion in the inner flux 
tubes of the model cloud. For At ^ 10 3 yr, ambipolar diffusion has become rapid 
enough to halt the growth of magnetic flux within this cell, (b) Same as (a), but for 
cell 5. The evolution is similar to that of cell 1 except that it exhibits a time delay, 
indicating that the region of piled-up magnetic flux is advancing outward from the 
center. 

Fig. 4.— Physical quantities in computational mesh cells I = 2,5, 10, 15,20,25, and 27, as 
functions of At. (a) Vertical (z) component of the magnetic field in the equatorial 
plane of the disk, normalized to -B Zj eq,co- Prior to the passage of the hydromagnetic 
disturbance (HMD), the magnetic field strength decreases because of inward advection 
of magnetic flux. At later times the magnetic field is increased by the compression 
induced by the outward-propagating HMD. (b) Column density, normalized to cr n)C o. 
In cells 2 and 5 the inflow of the neutrals is unaffected by the advancement of the 
HMD, and the column density decreases with time due to advection of matter into the 
central sink. In cells I > 8 the combination of the enhanced field strength (due to the 
piled-up magnetic flux) and the increased collisional coupling between the ions and 
neutrals at larger radii eventually leads to a deceleration of the infalling neutrals and 
the formation of a hydromagnetic shock. The column density increases immediately 
after the passage of the shock front; at later times the neutrals diffuse through the 
shock and are reaccelerated to free fall, with a n oc (At) s A t) —0.76 ^ s At ^ —0.55. 

Fig. 5.— Outward motion of the hydromagnetic disturbance as a function of At. (a) 
Instantaneous position of the HMD, normalized to the cloud radius R . (b) Speed 
of the disturbance (in the protostellar rest frame), normalized to C. Early on, the 
expansion of the HMD is nonsteady. At later times, the speed of the disturbance is 
subsonic, (c) Mach number relative to the infalling neutrals just upstream of the 
disturbance. Note that the time intervals between successive data points in this figure 
are larger than those in Figs. 2-4, reflecting a time-averaging of the motion of the 
HMD. For this reason the steady oscillations seen in Figs. 2-4 at later times do not 
appear in this figure. 

Fig. 6.— Spatial profiles of physical quantities in the typical model during the post-PMF 
epoch as functions ofr/Ro. Seven different times are plotted: At = 0, 2.17 x 10 2 yr, 
6.11 x 10 2 yr, 1.54 x 10 3 yr, 3.81 x 10 3 yr, 2.38 x 10 4 yr, and 1.48 x 10 5 yr. An 
asterisk on a curve locates, as in Fig. 1, the instantaneous position of the critical 
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magnetic flux tube, (a) Mass, in M . The protostar first forms at At m 0. (b) Vertical 
component of the magnetic field in the equatorial plane, normalized to -B 2j eq, c o- The 
dashed curve refers to a model in which the magnetic flux was assumed to be frozen 
into the neutrals at all times after PMF, shown at the last plotted time (Ate). The 
field strength is strongly amplified behind the expanding HMD for At > At2- For 
At 5 <; At <^ At 6 the field strength at t/Rq <; 2.5 x 1CT 4 is significantly reduced, 
reflecting the refreezing of magnetic flux into the infalling neutrals brought about by 
the greatly enhanced collisional coupling that occurs when v D > u Dcrit . (c) Neutral 
column density, normalized to <7 n , c o- The shock in the neutrals manifests itself by the 
enhancement and local maximum in the column density for At > At 4 . Far behind 
the shock the neutrals establish free-fall collapse, and o~ n oc r" 1 / 2 for r <^ 4 x 10 _5 i?o- 
The dashed curve refers to the model in which magnetic flux was frozen into the 
neutrals, at the same time as in (b). (d) Radial component of the magnetic field at 
the surface of the cloud, normalized to -B^eq.co- For At > At2, ambipolar diffusion 
halts the inward advection of flux and further bending of the field lines. The dashed 
curve refers to the frozen- flux model, at the same time as in (b). Continual inward 
advection of magnetic flux in the frozen-in model results in much greater bending 
of field lines than in the ambipolar-diffusion case, (e) Infall speed of the neutrals, 
normalized to C. The location of the hydromagnetic shock is marked by the abrupt 
transition from supersonic to subsonic inflow for At > At 4 . Free-fall collapse is 
established far behind the shock, and v n oc r -1 / 2 for r ^ 4 x 10~ 5 i?o — 35.5 AU. 
(f) Ion-neutral drift speed, normalized as in (e). In the region of rapid ambipolar 
diffusion behind the expanding hydromagnetic disturbance, v D ~ |i> n |. (g) Accretion 
rate, in M Myr" 1 . For At ^ At 4 the accretion rate is decreased behind the shock 
front, (h) Ratio of neutral-ion collisional timescale and gravitational contraction 
timescale. Near the inner boundary this ratio is ^> 1, indicating efficient ambipolar 
diffusion. Further away from the center, the degree of ionization increases and the 
magnitude of the gravitational field decreases, resulting in r ni /r gr < 1 and, therefore, a 
much better collisional coupling between the ion and neutral fluids, (i) Ratio of cloud 
vertical half-thickness Z(r) and radius r. By the time At 6 the tidal gravitational field 
of the central protostar has compressed the inner core to a disk that is geometrically 
thin. [Magnetic squeezing of the disk, not included in this calculation, further reduces 
Z{r).\ 

Fig. 7.— Mass of the gas surrounding the central protostar in the typical model, in M Q , as 
a function of r/R . All labels and times Atj are the same as in Fig. 6. 



Fig. 8.— Spatial profiles of quantities related to the stability of the core of the typical model 
cloud with respect to magnetic interchange, as a function ofr/R . All labels and times 
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Atj are the same as in Fig. 6. (a) Local mass-to- flux ratio, in units of the critical 
value for collapse. Unloading of mass from field lines due to rapid ambipolar diffusion 
behind the hydromagnetic disturbance results in positive values of d (In a n / B z ^ eq ) / dr 
inside the core at times At > At\. (b) Product of the linear instability growth rate 
and the neutral-ion collision time at At = At6 for the region of the core susceptible 
to magnetic interchange. Collisional coupling of the ions and neutrals is sufficient 
to allow the instability to grow in this region of the core, (c) Product of the linear 
instability growth rate and the kinematical timescale (= r/\v n \) for the same time 
and region of the core as in (b). An unstable mode can grow before being advected 
downstream by the fluid. 

Fig. 9.— Density profile of the typical model, normalized to n QtC0 , as a function ofr/R . All 
labels and times At,- are the same as in Fig. 6. In the inner core, n D oc r~ 2 throughout 
the post-PMF epoch. 

Fig. 10.— Physical quantities in the inner flux tubes of the frozen-flux model as functions 
ofr/R at seven different times At,-, (a) Ratio of ion-neutral drift speed and the 
neutral infall speed. In this model, V[ is set equal to v Q in the numerical code, and the 
drift speed v D is calculated from the ion force equation. As the core evolves, v D /\v n \ 
becomes ^> 1, indicating that the freezing of magnetic flux in the neutrals would 
break down, (b) Ratio of the r and z components of the magnetic field, (c) Ratio 
of the total radial magnetic force and the magnetic pressure force. After PMF, the 
ratios B r Z /B z eq and -F m ag,r/-^mag, P res scale as r~ 3 / 2 . The enhanced magnetic tension 
force following PMF would act to revitalize ambipolar diffusion during this phase of 
the collapse, which is indeed what is found to occur in the typical model presented in 
S 3.3. 
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